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PREFACE 


During the first week of April 1989, a workshop, entitled "Constitutive Relation- 
ships and Models in Continuum Theories of Multiphase Flows," was held at NASA's 
Marshall Space Flight Center. The purpose of this workshop was to open a dia- 
logue on the topic of constitutive relationships for the partial or per phase 
stresses, including the concept of solid phase "pressure" and the models used for 
the exchange of mass, momentum, and energy between the phases in a multiphase 
flow. This volume is the result of the stated objective of the workshop, namely: 
to provide a state-of-the-art report on constitutive relationships and models in 
continuum theories of multiphase flows. This written record is intended to be 
something of a "market place" for those engineers and scientists who are attempt- 
ing to implement multiphase flow theories for practical calculations. The focus 
on continuum theories was quite intentional, given the high level of engineering 
application that the continuum or "two-fluid" theory is seeing. Arguments per- 
taining to the applicability of the continuum approach versus other multiphase 
flow schemes, such as the Lagrangian-Eulerian or "tracking" techniques, are not 
broached. However, it is cautioned that even for certain broad categories of 
multiphase flows, the applicability of one class of theories over another is de- 
bated. Caveat emptor! 

The authors are from a variety of backgrounds and have, collectively, a knowledge 
of the spectrum of typical multiphase flows; i.e., from the dilute to the concen- 
trated in terms of the fraction or concentration of the dispersed phase. 

Within the general class of continuum theories, two (perhaps underexplored) 
formulations for multiphase flows are presented. These are the drift flux 
approach and the potential flow theory. These formulations provide a very useful 
and alternative perspective to the problem. However, as with the more familiar 
mass, momentum (with dissipation), and energy flux formulations, the drift flux 
and potential flow theories suffer from the problems of "closure" modeling re- 
quirements and validating the adjustable parameters. 

One author wrote in a summary of this workshop: "We are responsible to agree on 
the form of the equations, the averaging techniques to be used and the interpreta- 
tion of the terms in the averaged equations. Without this, we don't even have a 
common language to discuss the problem." The bulk of the presentations and 
discussion which resulted during the course of this workshop addressed this ob- 
jective. 

It is recognized from a physical perspective that the interaction of the dispersed 
phase particles (deformable or otherwise), amongst themselves, either via col- 
lision, "lubricated" near collisions, intermingling wakes or wake "drafting," are 
important processes. Hence, it would be an oversimplification to consider only 
single particle hydrodynamics when formulating constitutive relationships and 
momentum/energy exchange models. This requires modeling of paired and multi- 
body interactions with probability distribution functions (PDFs) and the solution of 
various integrals of products of PDFs over specified field spaces. This is an area 
where "borrowing" from the evolutionary structure of kinetic gas theory has 
proven to be useful. 

The analysis of dispersed phase stresses, specifically, the concept of solid phase 
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"pressure," is maturing. There are at least three mechanisms for normal or 
spherical momentum flux by a rigid dispersed phase: collisional transport, momen- 
tum transport due to the kinetic or deviation-about-the-mean motion of the 
particles, and the transport of hydrodynamic stresses which are present upon the 
skin of the particle. The first two mechanisms can be recognized as characteris- 
tics of the dispersed phase "cloud" or more specifically, as continuous, advectable 
fields of the dispersed phase. Structures, albeit complex ones, for the the formu- 
lation of these stresses are available and are presented in these proceedings. The 
latter mechanism is of unavoidable importance and requires analyzing the hydro- 
dynamical stress around a single particle and averaging these to achieve a con- 
tinuous field. Generally, the result is a dissipative (nonrecoverable) spherical 
stress, ascribed to the dispersed phase, containing, amongst other things, the slip 
velocity between the dispersed and continuous phases and the continuous phase 
viscosity. It must be cautioned that it is possible to include these same stresses 
redundantly in a model for the momentum exchange between the phases. 

Within this community, there is a growing consensus that the spin field, not to be 
confused with the vorticity field, of the dispersed phase is of importance. To 
date, attempts at modeling higher order interactions between the phases, such as 
lift forces, requiring some knowledge of the particle spin, have set that value 
equal to one half the local, continuous phase vorticity. This result is from single 
particle hydrodynamical analysis in an unbounded, linear shear flow and does not 
allow for any independent spin motion of the particles due to wall interactions, 
collisions, or "near" collisions. In addition, it couples the spin inertia of the par- 
ticles directly to the local vorticity field of the continuous field and does not 
allow for particle "spin-up" or "spin-down" in the presence of a local shearing 
flow. 

Finally, a suite of experiments which could be used for testing and validating/re- 
jecting various aspects of multiphase flow theory is sorely needed. In a manner 
similar to the single phase flow community, with their cavity flows and backward 
facing steps, this community needs to arrive at a consensus for a suite of experi- 
ments, so that when we finally all speak the same language, we can all shoot at 
the same target. 
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ABSTRACT 


The main characteristics and the potential advantages of 
generalized drift flux models are recalled. In particular it is 
stressed that the issue on the propagation properties and on the 
mathematical nature (hyperbolic or not) of the model and the 
problem of closure are easier to tackle than in two-fluid models. 

The problem of identifying the differential void-drift closure law 
inherent to generalized drift flux models is then addressed. Such 
a void-drift closure, based on wave properties, is proposed for 

bubbly flows. It involves a drift relaxation time which is of the 
order of 0.25 s. 


It is observed that, although wave properties provide essential 
closure validity tests, they do not represent an easily usable 
source of quantitative information on the closure laws. 
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I. INTRODUCTION 


Generalized drift flux models were recently shown (Boure, 1988b) 
to be attractive alternatives for current 1-D two. fluid models. 
Drift flux models are characterized by the uses of a single 
momentum balance (the mixture balance instead of two phasic 
balances) and of a void-drift closure law. In classical drift flux 
models the void . drift closure law is expressed through an 
algebraic equation, which amounts to ignoring nonequilibrium drift 
effects. In generalized drift flux models, the void-drift closure 
equation is a partial differential equation. 

Generalized drift flux models and two- fluid models are compared in 
the next section and in table 1. The comparison brings out the 
drawbacks of two-fluid models. However, both kinds of model must 
be complemented by closure laws. In particular, generalized drift 
flux models need a void . drift closure law which remains to be 
specified . 

Since generalized drift flux models were introduced to account for 
the properties of kinematic waves (Boure, 1988a), it seems 
logical to use these properties to identify (i.e. to evaluate the 
coefficients of) the void-drift closure equation. The purpose of 
the present paper is to discuss the identification problem, 
assuming that the void- drift closure equation may be approached by 
a quasi linear differential equation of the first order. 

!X. A REMINDER ON GENERALIZE D DRIFT FLUX — MODELS — AND TWO -FLU ID , 
MODELS 


The comparison between generalized drift flux models and two- fluid 
models is summarized in table 1 which, like must of the substance 
of this section, is taken from Boure (1988b). The mass and energy 
balances are parts of both kinds of models, and they are not 
discussed further hereafter. It is only noted for completeness 
that they require a few closure equations, in particular for the 
mass and energy transfers at the walls and at the interfaces. 
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The momentum balances contain the two phasic averaged pressures 
and P L , the subscripts G and L corresponding to the two phases. 
For the following discussion, it is convenient to express P Q and 
P L in terms of the average pressure P and the pressure difference 
P LG , with : 


P £ a P G + (1 - a) P L , P LG = P G - P L (2.1) 

a being the void fraction. Now, in the set of balance equations, 
the two phasic momentum balances are equivalent to a subset of two 
equations, namely the mixture momentum balance and the pressure 
difference equation , obtained on eliminating P between the two 
phasic balances. 


The mixture momentum balance may be written : 

dP ( d W G d w g^ fi W d W. 

_ + c Po ♦ w 0 _ j ♦ (1 - a) p L (_ + „ L ^-i 

( 2 . 2 ) 


3R dr 

+ dl + dl + *^ (W G * W l) - la + F w + PS = 0 


the terms of the second line representing respectively : 

a term accounting for fluctuation and transverse distribution 
effects 

a term accounting for longitudinal stress variations 

- an interfacial mass transfer term 

- a surface tension term (interfacial) 

- a friction term (walls) 

- a gravity term. 


t and z are respectively the time and the space variables o, and 
P L are the phase densities, W Q and W L are the phase average 
velocities, g is the gravity acceleration, and : 

P = “ P G + (1 - a) P L (2.3) 
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w being the local instantaneous velocity along Oz and j the local 

k ... 

instantaneous deviatoric stress tensor, indicating the 

conditional time or ensemble averaging operator and < > the 

space averaging operator, and n being the unit vector of the Oz 

z 

axis, R and r are defined as : 


G L 

R 5 (oc p G (w G - W Q ) 2 + (1 - a) p L (w L - W L ) 2 y (2.4) 

r = ~ ^ |ar G + (l-a)x L j*D^*I} (2.5) 

Independently of the mass transfer, already present in the mass 
balances, eq. 2.2 requires four closure laws for R, t, I ff and F y . 
Turbulence effects are present through R. 


The pressure difference equation may be written : 

3 W r \ /d W, 


a (1 - a) 


[ 


fi «0 0 "o'! 

' ?T + W » — J " 


/d w L a w L \ 

+ W. — — 

{ dt 1 dz ) 


3 P LG 

+ 


3z 


+ Lr + L r + M lg [(1 - a) (W G - W GI ) + a (W L - W LI )] 


( 2 . 6 ) 


+ L,j + L F y + Lp j a (1 oc) p G ^ g 0 


the terms of the second and third lines representing respectively 

- a term accounting for fluctuation and transverse distribution 
effects 

- a term accounting for longitudinal stress variations 

- an interfacial mass transfer term (W QI and W L j are interfacial 
averages of the phasic velocities) 

- an induced inertia ("added mass") term 

- two friction terms (respectively wall and interfaces) 

- a gravity term, with 


Pgl - Pl “ Pg 


(2.7) 


The above terms are defined in Boure (1988b) . For instance 


o 


( 2 . 8 ) 


** = U-a) ^ ( a P G ( W G - W G ) 2 ) * a ^ (d-a) P L ( W L - w l > 2 ) 

L t a - (1-0.) 2L(« f‘- n ; ) • n ; a a 1. (d-«) l' • Sj ) • ^ (2.9) 

Equation (2.6) requires seven closure laws for L_ L W 

I' * T” F G I 1 

W L I ’ L i » l fw» i * Moreover, retaining eq. (2.6) implies, at least 
in theory, retaining also the interfacial momentum balance, which 
requires one supplementary closure law. 

In two-fluid models, eqs (2.2) and (2.6) are both implicitly 
written. A crucial point is that eq . (2.6) is stiff. This means 
that, due to the relative orders of magnitude of its terms, small 
variations on the closure laws (and in particular on the induced 
inertia and interfacial friction closure laws) induce large 
variations on P LG (when P L Q is not merely ignored) and/or W Q . 
Hence the difficulty to adjust the very badly known closure laws 
of eq. (2.6) to avoid unrealistic values of the light phase 
velocity . 

A second crucial point is that, whenever eq. (2.6) is closed by 
algebraic laws only, a non-hyperbolic set results : two 
characteristic velocities are complex conjugate, with the 
consequence that the model is unconditionally unstable. 

In drift flux models, eq. (2.6) is not written, which is 
acceptable since the value of P^ G does not really matter in 
practice. The two foregoing difficulties are not encountered. 

Besides the closure laws already mentioned, current two- fluid 
models are closed through an assumption on P L G (the set of 
6 phasic balance equations involves 7 dependent variables, namely 
two pressures, two velocities, two enthalpies and the void 
fraction). Such an assumption imposes an artificial constraint on 
the pressures and pressure gradients. It disturbs the description 
of the corresponding propagation phenomena. 
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Generalized drift flux models are closed through the direct 
description of the void-drift dynamic dependency. Assuming a 
quasi -linear , partial differential relationship of the first 
order, and using the convenient variables : 

W ^ a W G + (1 - a) W L , 6 £ a <1 - a) (W G - W L ) = a (W G - W) (2.10) 


(center of volume velocity and drift), it may be approached by . 




doc 

~ m * 


as 

+ — + 

at 


w. 


as 

dz 



( 2 . 11 ) 


Equation (2.11) requires seven closure laws, i.e. the same number 
as eq. (2.6) but with simpler physical significances and more 
straightforward consequences : f is the fully-developed drift 
value, 6 a relaxation time, I and II are respectively the sum and 
product of the two characteristic velocities corresponding to the 
kinematic waves, t, expresses inertia effects, W 2 and W 4 are 
averaged velocities close to W. 

It can be concluded that developing a closure set for eq . (2.11) 
and using generalized drift flux models appear as less hazardous 
and hopefully easier than developing a closure set for eq. (2.6) 
and using current two-fluid models. 


m . TWVT.ITK.NGE OF THE VOID - D R TFT CLOSURE ON THE PROPERTIES O F 
KINEMATIC WAVES ("Direct problem") 


As long as the kinematic wave velocities are small with respect to 
the sonic velocity, the properties of small harmonic kinematic 

waves of the form 


7 


X 


i (cot - kz) 


(3.1) 


= x o e 

where x Q (constant) CO and k are real or complex quantities, result 
from the approximate dispersion equation (Boure, 1988b) 

<0 - C a k + i 0 (co 2 - I cok + II k 2 ) = 0 (3.2) 


where : 


C a = W + f a# 



(3.3) 


A first consequence is that the kinematic wave properties, which 
do not significantly depend on tj, W 2 , W 4 , cannot be used to 
evaluate these quantities (r| is related to the sonic velocity, W 2 
and W 4 have only weak influences) . 

In the exploitable experiments (Tournaire, 1987, Boure, 1988a) co 
is imposed and the kinematic wave velocities V and their spatial 
amplification coefficients k* result from the data processing. 
Eqs . (3.1) and (3.2) must therefore be used with CO real and : 

k = k f + i kj (3.4) 


from which 


k: z i (cot - k r z) 

x = x Q e e r 


(3.5) 


V 



(3.6) 


k p , which does not depend on the frame of reference, can be used 
instead of co, which does, to characterize a wave. 


Introducing eqs. (3.4) and (3.6) in the dispersion equation (3.2) 
and separating the real and imaginary parts lead to : 
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k r [V - c a + 0 Z V kj - 2 0 n kj ] = 0 


(3.7) 


- c a k, + 0 V 2 kj; - 0 Z V kj! + 0 n (kj: - k?) = 0 (3.8) 

Equations (3.7) and (3.8) provide the solutions to the direct 
problem, viz. computing the properties of the kinematic waves of 
wave number k p when C a , 0, X and IT are known. There are two 
solutions corresponding to two modes (noted with the subscripts 3 
and 4) . In the experiments it was found that modes 3 and 4 are 
respectively predominant at low void fractions (0 < oc < 0.25) and 
at "large" void fractions (a > 0.30) 


In particular for k p =0, the two solutions are : 


k i 6 < C a - C 3> ( C oc - C 4> 
k _ q with — - 


V = C 


(3.9) 


and 


kj = 


0 n ' v c<x 


1 ^ 1 1 1 

C ; C 3 C 4 C <x 


(3.10) 


C 3 and being defined by : 


2 = C 3 + C 4 


n s c 3 c 4 


(3.11) 


For k f oo the two solutions are : 


V 


C 


a 


C 3 


0 C 3 (C 4 - C 3 ) 


V = C, 


k i = " 


C 4 - c a 


0 c 4 (C 4 - C 3 ) 


(3.12) 


(3.13) 
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In equations (3.7) and (3.8) the three quantities 9 k p , k- are 
present only through the two products (0 k p ) and (9 k- ) . In actual 
experimental runs, k p is never zero and the equations may be 
written : 

v - c a + ( c 3 + c 4 ) v (9 k ,. ) - 2 c 3 c 4 (6 kj ) = 0 (3.14) 

9 k, (0 k ( ) 2 

- C a + V 2 - (C 3 + C 4 ) V + C 3 C 4 - C 3 C 4 = 0 (3.15) 

(0 k r ) 2 (0 k r ) 2 

Since in the exploitable experimental data, V does not 
significantly depend on k p (no^ significant dispersion) , it is 
convenient to eliminate V between eqs . (3.14) and (3.15). Eq . 
( 3 . 14) yields : 


C a + 2 C 3 C 4 (9 k f ) 

V - (3.16) 

1 + <c 3 + C 4 ) (9 k- ) 


and eq. (3.15) yields : 


V = 


C 3 + C 4 


(C, - C,) 


+ C„ 


0 kj (0 kj ) 

+ Ci c. 


(0 k r )^ 


4 


(3.17) 


(0 k p ) f 


The solutions for 0 k- are then the real solutions of the 
equation 


2 C a - (C 3 + C 4 ) - (C 4 - C 3 ) (0 kj ) 

1 + (C 3 + C 4 ) 0 kj 


(C 4 - C 3 ) c + 4 C c 


0 kj (0 kj y 

+ 4 C, C, 


(0 k p )‘ 


J 3 4 


(0 k r y 


(3.18) 


Computing directly 0 kj as a function of 0 k p from eq. (3.18) is 
not straightforward. It is more convenient to transform eq. (3.18) 
to express 0 k f as a function of 0 kj 
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(3.19) 


(0 k r ) 2 

0 k i 


[c* + c 3 c 4 e k,. ] [l + (c 3 + c 4 ) e k, ] 2 

( C a“ C 3) (Coe" C 4> " < C 4- C 3) 20 k i l C a + C 3 C 4 0 k i 1 


For a given value of k- , eq. (3.19) yields zero or one real 
positive value of k p . 


IV. IDENTIFICATION OF THE VOID-DRIFT CLOSURE EQUATION FROM 
KINEMATIC WAVE PROPERTIES ("Inverse problem") 


The problem posed in this paper is the determination of C a , 9, Z 
and IT, using the experimental data on kinematic waves. It is the 
inverse of the problem of section 3. 


The principle of the method is to write eqs . (3.7) and (3.8), for 
instance, for several sets of experimental conditions for which 
k p , V and k j are known and to use the resulting equations to 
compute C a , 9, Z and IT. 

In the exploitable data, as already noted, V does not depend 
significantly on k p . Accordingly, when a single mode is 
predominant, V may be expected to be close to both C-j (or ) and 
C a (or C^) . This has two consequences : 

1. Whenever a single mode is predominant, the experimental 
correlation for V should be a correlation for C a as well. This is 
corroborated by the fact that is satisfies the definition (3.3). 
In the experimental conditions of Tournaire (1987) (upward 
vertical flow, low pressure), it leads to (Boure, 1988a, f and 
C a - W in m/s) : 


For 0 < a < 0.2 (mode 3 predominant) 

f = 0.22 a (1 - a) [1 - 1.25 a (1 - a)] 

C^ - W = 0.22 (1 - 2a) [1 - 2.5 a (1 - a)] 


(4.1) 
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For 0.3 < a < 0.41 (mode 4 predominant) 


f = 0.22 a. - 0.028 

C a - W = 0.22 


(4.2) 


For 0.2 < a < 0.3 (the two modes coexist) the values of f and 
C a - W may be interpolated between (4.1) and (4.2) with, from the 
experimental data : 


C a - W a 0.08 (4.3) 

for a = 0.25 

2. When mode 3 (respectively mode 4) is predominant, C a - C 3 
(respectively C 4 - C a ) should be "small". 

being known and V eliminated, the problem may now be 
reformulated, eqs . (3.7) and (3.8) being replaced by eq . (3.18) or 
(3.19) to be solved for 9, C 3 , C 4 . Only mode 3 results are 
exploitable since mode 4 results for k- are too few in number. In 
view of the forms of eqs. (3.18) and (3.19), the foregoing problem 
is far from simple. This is confirmed by fig. 1 in which, as 
su ggested by eq. (3.19), the experimental results for - kj;/kj are 
plotted as a function of — kj , and which exhibits an important 
scatter (in the representation of fig. 1, the points corresponding 
to Ikj I < 0.1 are subject to large errors and therefore 
meaningless. They were not plotted in the figure). 

In the domain in which mode 3 is predominant, the conditions : 

C 3 < C oc < C a < C 4 (4-4) 

may be expected to hold (see Boure , 1988a). They entail 

C a ~ C 3 ; 1 ^ C 4 - C oc Coc 

0 C 3 ( C 4 — C 3 ) 9 ( Cj + C 4 ) 9 C 4 (C 4 ~ C 3 ) 9 C 4 
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Then, as k ; decreases from zero, kji/kj resulting from eq. (3.19) 
starts from the value corresponding to eq. (3.9) with a slope 
which may be of either sign. However, when Ik^ I is sufficient, 
k^/kj also decreases. It tends towards - oo for the value 
corresponding to eq. (3.12). For the values of k- comprised 
between those given by eqs . (3.12) and (3.13), there is no 
physical solution (kj: < o) . Finally, for the values of k ( 
comprised between those given by eqs. (3.13) and (3.10), there is 
a solution again corresponding to mode 4. 

By trial and error, 9, C 3 and C 4 may be adjusted to fit the curve 
representing eq. (3.19) to mode 3 data. In view of the absence of 
dispersion C 3 may be expected to be close to C a . 

For 0 < a < 0.20 : the following set of values is acceptable 

9 = 0.25 s ) 

C a - c 3 =0.02 m/s > (4.6) 
C 4 - C a = 0.08 m/sJ 

Precise adjustement would need more accurate data for the wave 
velocities and especially for the damping/amplif ication 
coefficients. Such data does not exist and cannot be expected to 
be obtained soon in view of the available instrumentation. For 
a > 0.25, no sufficiently accurate data is available. 

Equations (4.1) and (4.6) confirm that, even at low pressure, the 
relevant velocity differences corresponding to fully developed 
conditions are fairly small in bubbly flows. They are negligible 
as soon as W is large enough (say 3 m/s) . Accordingly the 
correlations for C a - W, C 3 - W, C 4 - W are probably not crucial. 
On the other hand, the drift relaxation time 9 is an essential 
parameter of the generalized drift flux model. 
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V. CONCLUSIONS 


After a reminder on the main characteristics and the potential 
advantages of generalized drift flux models , the problem of the 
identification of the differential void-drift closure law they 
imply has been addressed. 

Two advantages of generalized drift flux models versus complete 
two- fluid models are : 

1. The correct description of the kinematic wave phenomena and the 
straightforward control and interpretation of the mathematical 
nature (hyperbolic or not) of the model set of partial 
differential equations that they enable. 

2. The simplification of the closure problem, involving only 
closure laws of simple physical significance and easy to assess. 

On the other hand, it has been found that the available data on 
kinematic wave properties is not quite adequate to enable the 
identification of the void-drift closure law. Such an 
identification would need a very good accuracy (difficult to reach 
in practice) on the wave damping or amplification coefficients. 
Accordingly kinematic wave properties seem to be more useful as a 
closure validity test than as a source of quantitative information 
on the closure laws. 

However, a void- drift closure, suitable for bubbly flows, has been 
adjusted on the available kinematic wave data. It involves several 
velocities which differ only slightly from each other and from the 
average fluid velocity but whichare necessary to the description 
of the kinematic wave properties. It also involves a drift 
relaxation time which is an essential parameter of generalized 
drift flux models and which was tentatively found to be of the 
order of 0.25 s for the upward flow of air-water mixtures at low 
pressure . 
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Table 1 

Simplified comparison of modeling strategies 


Drift flux approach 


Common features 


Two-fluid approach 
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NOMENCLATURE 


d£, ds 

E 

f 

g 

I 

j 

k 

m 

m 

M? 

P 

P 

P 

R 

t 

u 

v 

w 

a 


V 

7T 

P 

$ 


Elements of length, area 

Exertia, (cti/3 — 1) 

Force per unit volume of a phase 
Body force per unit mass 
Unit dyadic 
Volumetric flux density 
Kinetic energy density 
Geurst’s added mass coefficient 
Momentum density 
Added mass force 
Mean pressure in a phase 
Macroscopic pressure 

Combined momentum flux and stress tensor 

Density ratio pi/p 2 

Time 

Microscopic velocity 
Average velocity 
Relative velocity Vj — V 2 
Volumetric fraction of a phase 
“Resistivity,” (1) 

Macroscopic potential 
Specific momentum, (24), (25) 

Density 

Macroscopic potential 
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Subscripts 


With particles at rest 
Phase 1 
Phase 2 
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INTRODUCTION 


Potential flow theory for a single fluid has been established for many years. Although 
its limitations for describing real motions axe well known, it does provide a self-consistent 
structure for analysis and often provides a reasonably accurate description of at least part 
of the flow field of fluids with low viscosity. 

As an example of two-phase flow, one can imagine a suspension of particles in a fluid 
that obeys all the requirements of potential theory at the microscopic level If the only 
forces acting on the particles are conservative,” it would appear that their motion might 
also reasonably be expected to be describable in terms of a suitable potential Averaging 
of these potentials would lead to macroscopic potentials, true properties of the mixture, 
that should be related in some way to the average motion of the phases. Indeed, previous 
attempts to determine the inertial coupling terms in the two-fluid model have implicitly 
assumed potential flow at the microscopic level 

A complete two-phase theory of this type will be as idealized as was the classical 
theory of single-phase potential flow. However, it should be useful in the same sorts of 
ways, both as an approximation in many situations and as a standard that must at some 
level be consistent with other approaches, such as those that attempt to define closure 
relations for a set of averaged basic equations. 

This paper describes some features of two recent approaches along these lines. The first 
is based on a set of progressive examples that can be analyzed using common techniques, 
such as conservation laws, and taken together appear to lead in the direction of a general 
theory, the tactic used in [1], The second is based on variational methods , a classical 
approach to conservative mechanical systems that has a respectable history of application 
to single phase flows. The latter approach, exemplified by several recent papers by Geurst 
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[2-5], appears generally to be consistent with the former, at least in those cases for which 
it has been possible to obtain comparable results. 

PROGRESSIVE EXAMPLES 

The theory developed in [1] starts with a situation where the particles are at rest. 
Fluid flows past these particles as though a porous medium . At the microscopic level 
the equipotentials are not “smooth” but they do not differ much from the more gross 
macroscopic equipotential. Differences between the two levels of equipotential are smeared 
out over lengths comparable with the particle size. I do not have a rigorous proof, but I 
believe it is valid to treat these two equipotential surfaces as essentially identical for most 
purposes. This is not true of the other properties, such as pressure and velocity, that vary 
more at the microscopic than the macroscopic scale and must be averaged carefully. 

Just as in electrical conduction past a matrix of non-conducting particles, the macro- 
scopic fluid flux will be proportional to the macroscopic potential gradient, the “resistivity” 
being represented by a factor j3 that depends only on the particle geometry and the void 
fraction, as long as the arrangement is isotropic. We therefore have 

jo = (1) 

Since f3 is unity for unimpeded flow, it will be greater than one when particles are 
present. A crude description of the situation could be to say that some of the fluid is held 
up” or “entrained” by the particles so that only a fraction of the space is available for 

direct flow. 

The average velocity of the fluid is the relative velocity , 

w = ^ = (2) 

<*i <*i P 
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The kinetic energy of the fluid in unit volume confined between two equipotentials is, 
from a standard theorem of potential flow. 


— ~ 2 Pl J <f> u ' ds = — -pi -j 0 = -p\ot\ftw 2 


(3) 


The kinetic energy per unit volume of phase 1 is 


1 *0 1 9 / n\ 

z 


(4) 


If we denote the fluid velocity at the microscopic level by Ui , (2) and (4) are equivalent 


to 


w =< Ul > (5) 

w 2 ai/3=< Ui > (6) 

Invoking the Schwarz Inequality, it is clear from (5) and (6) that ct\fi is greater than 1, a 
proof pointed out to me by my student Chao Luo. 

By considering the changes in the fluid kinetic energy resulting from a uniform volume 
change for every particle it can be shown [1] that the difference is mean pressure between 
the phases is 


1 2 d , „ x 

P2 - Pi = 2? lW Ql - 1) 


(7) 


The additional “1” in (7) appears to be gratuitous. It is introduced at this stage 
because the “exertia,” defined [1] as 
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E = (otifi — 1 ) 


( 8 ) 


turns out to reappear in numerous contexts and to be one version of an added mass 
coefficient. 

If we now superimpose a uniform velocity V 2 on the above motion, this is equivalent 
to superimposing a corresponding potential gradient and we obtain 


Vj + £w -= — 


( 9 ) 


where the relative velocity is 


w = Vi — V2 


( 10 ) 


The net momentum and kinetic energy densities are then 


m = + P2OL2V2 


(ii) 


k = 


-piaiUj + ~P20i 2 vl + -pia x Ew^ 


( 12 ) 


From (11) and (12) we may deduce that the effective equations of motion of a uniform 
suspension accelerating under the influence of a uniform macroscopic pressure gradient 
and body force fields are 
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V] + Ew = 


VP 


ai xw 


(13) 


Pi 


+ gi + 


Pi 


Pl^l ^ VP , _ , a 2 XW 

V 2 E W — h g2 H 

P2<*2 P2 P2 


aj and a 2 are arbitrary vectors subject to the constraint 


(14) 


(oiai + o 2 a 2 )xw = 0 (15) 

It would appear that the assembly must be incompressible if it is to be truly “uniform” 
in a pressure gradient. If gi is a conservative force field, comparison between (9) and (13) 
would seem to indicate that Vxai x w is zero and therefore Vxa 2 x w is as well from (15). 
Eq. (14) then suggests that if g 2 is conservative, the combination on the left-hand side is 
the gradient of another “potential” which we could define in a form somewhat resembling 

( 9 ) : 


V2 _ = — V 77 (16) 

P 2<*2 

A different argument is used to derive the one- dimensional equivalent of (16) in [ 1 ], In each 
case, the interpretation of $ and rj could be, in the classical view, in terms of impulsive 
pressures and body forces necessary to set up the motion. This derivation, made for a 
uniform suspension, requires further argument, or a leap of faith, if it is to be applied 
more generally to a dispersion in which w and e*i, and hence /3 or E, vary with position. 

Eqs. (13) and (14) show that the exeriia plays the classical role of a coefficient of 
apparent mass , being proportional to various alternative definitions [l] in the literature. 
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Some other results derived from mechanistic arguments [1] axe: 

• The combined momentum flux and stress tensor for the suspension: 


p = <*i(pi V 1 V 1 + Pi I) + « 2 (P 2 V 2 V 2 + P 2 I) + aipi-Eww ( 17 ) 

• Bernoulli’s equation for fluid flowing steadily past a stationary particle matrix: 


Pi = Poi ~^Piw 2 (l + E) (18) 

• The force per unit volume on a particle in a stationary lattice: 

f 2 = v P2 = vQ ft »M^-) (19) 

When these results are used to check several hypothesized forms of the equations of 
motion in the two- fluid model of two-phase flow, there are found to be discrepancies [7], 
except when Geurst’s equations are used. 

GEURST’S EQUATIONS 

In a series of papers [2-6] Geurst has used variational methods to derive a number of 
results, starting from the hypothesis that the kinetic energy per unit total volume is 


k = 


\pi a i v \ + ^P20t 2 v 2 + ^Pimw 2 


which is the same as (12) with the alternative definition 


( 20 ) 


rn = — 1 ) = ot\E 


( 21 ) 
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Initially Geurst [2,3] developed a one- dimensional theory, introducing Lagrange mul- 
tipliers that are the same as the potentials $ and 7 / in (9) and (16). He went on to 
derive equations of motion, momentum, fluxes, pressures, etc. that are compatible with 
the results in Section 2, where there is a clear equivalence. 

In three-dimensions the derivation is complicated by the introduction of Clebsch po- 
tentials and Lin constraints. Details are not provided in earlier publications [2-4] but they 
do appear in a recent one [6] in which the equivalents of (9) and (16) are expressed sis 


TTi = + V’lVxi (22) 

7r 2 = V<^2 + ^2VX2 + — V^ n (23) 

P 2 

where TV\ and 7r 2 are generalized specific momenta defined as 

TYl 

ir 1= vi (v 2 -vj) (24) 

at 1 

^r 2 — v 2 -(- 2 Vj) (25) 

P 2«2 

which are identical with the left-hand sides of (9) and (16), in view of (21). 

The (j )’’ s and \’s are five “potentials” or Lagrange multipliers and n is the bubble 
density. The detailed interpretation of these terms is not important in the present context. 

Geurst’s equations of motion appear in many equivalent forms. The versions selected 
in [8] may be written in the present notation as 
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— (/SlttlXi) + V • (/?!«! V i^i) = — Ofi Vp2 - V 




-pi(m + aim )w 


+oifi + pim{x 2 - v x ) • (Vv 2 )^ 


( 26 ) 


— (/? 2 a 2 ^ 2 ) + v • (/>2Q'2V 2 7r2) = -a 2 Vp 2 + a 2 f 2 - Pim(v 2 - Vi) • (Vv 2 ) r ( 27 ) 

which have the form of conservation laws for 7Ti and k 2 . m' denotes dm /da 2 . 

The two continuity equations are 

d 

^(/>i«i) + v • (^xaiVx) = 0 (28) 

d 

-^(p 2 a 2 ) + V • (p 2 a 2 v 2 ) = 0 (29) 

If (7) is accepted and (21) is used, we find that the term involving w 2 on the right-hand 
side of (26) may be rewritten as 

]-p!(m + aim')w 2 = ai(pi - p 2 ) (30) 

Using (28) through (30) it is straightforward to arrange (26) and (27) into the forms 

^(irj) + V (vx -7Ti - (V 2 - Vx) 2 ^ -VxXVXTTx = -^! + — (31) 

dt \ 2 2«x ) pi pi 

2 ) 4- V (v 2 • 7T 2 - ^ - V 2 xVX7T2 = — — + — (32) 

dt V 2 / P 2 P 2 
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Eqs. (31) and (32) are versions of the Bernoulli Equations derived in [8] without 
identifying Tt\ and tt 2 with gradients of potentials. 

Now, if either pi is constant or Vpi is parallel to Vpi, and V x (fi/pi) = 0, we may 
take the curl of (31) and obtain 

p\ 

— (VXTTj) — VxV] x(VX7Ti) = 0 (33) 

at 

which is a conservation law for the vorticity of . If we consider the rate of change of the 
net vorticity of 7Ti threading a loop moving with the velocity Vj, we obtain 

= j £(Vx7T 1 ) ds + J (Vxnq) • (vi xd^) (34) 

The two terms in (34) represent the sum of in-place changes in the flux through the loop 
and contributions picked up by the motion of boundaries. Changing the order of the scalar 
triple product in the final term and invoking Stokes’ Theorem, we obtain a surface integral 
of (33) which is identically zero, in other words, 

1 /”'" = ° < 35) 
This result was obtained by Geurst [3]. A similar conclusion follows for (Vxtt 2 ). 

Now, if both phases come from a region in which there is no curl to TT\ and tt 2 (for 
example, a stagnation region), conservative body forces act and density gradients (if any) 
are parallel to pressure gradients, then throughout the flow, in view of (35), it should be 
true that 
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V x 7T] = 0; 


Vx«-2 = 0 


(36) 


which implies that both 7Ti and 7T2 are gradients of suitable potentials and only the first 
term is needed on the right-hand sides of (22) and (23) which reduce to the form of (9) and 
(16). This development apparently makes more explicit the conditions for the existence of 
two-phase potential flow , requirements which parallel those for the classical single-phase 
case. 

DISCUSSION 

The previous sections of this paper have outlined two approaches to potential two- 
phase flow. Each has a justifiable theoretical base and is self-consistent. Moreover, both 
approaches appear to give the “right” prediction for several well-defined situations [7] while 
some other formulations fail these tests. 

In order for these ideas to blossom further there need to be: 

a) Further developments, from the same basic set of assumptions, that encompass 
more generality. 

b) More rigorous derivations that clearly explain the order of approximation involved 
in treating the flow of discrete entities as a continuum. 

c) Reconciliation with alternative approaches, particularly those involving averag- 
ing. 

d) More solutions to specific problems that can be thoroughly investigated for con- 
sistency. 

e) An understanding of outstanding incompatibilities between these approaches and 
various other theories, with a clear explanation of what has “gone wrong.” 
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I hope that some of these items can be discussed at this meeting. At this time I only 
wish to briefly address (d) and (e). 

Apart from the “tests” described in [1] and [7] and perturbation techniques leading 
to the description of wave propagation and stability in bubbly flow [2,3,5], I know of 
only two complete solutions to Geurst’s equations. The first involves steady flow of two 
incompressible fluids from a common stagnation region. The simple conclusion that is 
reached is that the void fraction is constant and the velocities of both phases are the 
gradients of potentials that are proportional to each other and may be borrowed from an 
equivalent single-phase flow. The velocity ratio, or slip ratio is 


V2 

Vi 



(37) 


where 


R = p\! p2 


(38) 


In the above derivation use was made of Maxwell’s [9] approximate expression for the 
exertii x, 


E » y (39) 

The second solution again uses (39) and leads to a relationship between the velocities 
in unsteady incompressible flow [8] 
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This has the same form as a relationship which would be obtained by eliminating the 
pressures from equations of motion that ignored inertial coupling, except that the effective 
densities are changed. There is opportunity to test (40) by comparison with the one- 
dimensional transient response of fluidized beds. 

Regarding item (e) I will repeat, in a slightly different form, Geurst’s equations of 
motion [4] as presented in [8] 

— (piaivi) + V • (piaivivi + pimww) = 

dt 

— aiVpi - (pi ~P2)Vai + Mf + aifi 


— (P 2 « 2 v 2 ) + V • (p 2 a 2 v 2 v 2 ) = -a 2 Vp 2 - + a 2 f 2 

dt 

with the added mass term expressed as 


(42) 


Mi = p\m 


■(| + V2 .v)v 2 -(| + V ,v) 


Vi + w X V X Vi 


W‘ 


+/5imV — + w 


Q 

— (pim) + V • p\m\2 

at 


(43) 


This set of equations contains several more terms than one would find in most similar 
expressions in the literature. By dint of these extra terms, several “tests” are passed that 
other formulations fail [7). Specifically, these tests involve (18) and (19), which reduces 
for a dilute dispersion to Taylor’s expression [10) for the force on a stationary object in 
an accelerating flow. It would be desirable to devise other “tests” that might help to 
discriminate further between true and false expressions. 
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ABSTRACT 


In this study, a self-consistent derivation of the conservation laws is given for flows of 
a fluid-solid mixture. A unified analytical framework for obtaining constitutive relations is 
provided This analysis uses a control volume/control surface approach that is widely used 
Fn fluid mechanics. All terms in the governing equations and the constitutive relations are 
written in terms of the mass-weighted averages except solid concentration. It is believed 
strongly that the mass- weighted average is the natural bridge between 
constitutive relations. The derived momentum equations contain terms that differ from 
existing models except that of Prosperetti and Jones (1984). However, their assumptions 
are nof needed here. Special attention is given to the solid phase pressure. The physical 
KriSof previously assumed form for this pressure (Givler 1987) becomes clear. A number 
of related phenomena are also discussed. These include the anti-diffusion and anisotropic 
normal stresses. The energy equations are also different from existing models. But detail 
discussion on the energy equations is left to future work. 


I. INTRODUCTION 

Modeling a flowing fluid-solid mixture starts from writing down a set of governing 
equations. These equations describe the conservation of mass, momentum and energy. In 
the early stage, the popular approach was to view the mixture as a single phase mate- 
rial. Consequently, the following type of equations were used for mass and momentum 
conservations (see Zuber 1964, Ishii 1975, 1977). 

%^ + V(p m )u m = 0 (1) 

at 

S- (Pm) u m + ^ • (Pm) u m u m = V • T + (Pm)g (2) 

at 

In which, ( ) represents the ensemble average, p m is the mixture density, u m is the mixture 
velocity, T m is the mixture stress and g is the body force per unit mass. 
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Recently, the two-phase flow approach becomes more popular. Using this approach, 
the conservation equations axe formulated for each individual phase separately. This ap- 
proach allows for direct modeling of a fluid-solid mixture when the two phases have dis- 
tinctly different dynamics. For mass and momentum, these conservation equations are 
generically written in the following way (see Ishii 1975, Drew 1976, 1983, McTigue et al. 

1986). 

d(p 3 c)/dt + V • (p,cu) = 0 ^ 

&(PfO- - c ))/dt + V • (p/(l - c)v) = 0 (4) 

d(p a cu)/ dt + V • (p M cuu) = <p, C g) + (m) + V • (cT.) 

d (p f (l - c)v)/dt + V • (p f ( 1 - c)vv) = ( Pf ( l _ c ) g ) - ( m ) + V • ((1 - c)Tf) (6) 

where p a and p f are densities of the solid and fluid phases, u and v are the velocities of 
the solid and fluid phases, m is the interaction force per unit volume of the mixture, c is 
the local solid concentration (equals to 1 at a solid point and 0 at a fluid point), T. and 
T f are the solid and fluid stress respectively. The above equations do not consider the 
phase changes at the interface. The energy conservation equations have not been studied 
as extensively as the mass and momentum conservation equations. 

Because the two phases are separated, available information on a single particle’s mo- 
tion and the particle-particle interactions in a fluid environment are incorporated directly. 
In a mixture model, these informations will first be utilized to obtain the drift flux term in 
T m ° f Eq ‘ Preference of the two approaches apparently depends on the application. 
However, since the mixture model can be derived from adding together the two-phase 
equations, the two-phase approach is considered more fundamental. 

Many mathematical models have been derived based on Eqs. (3)-(6). There are two 
issues in modeling terms appeared in those conservation equations. First, what kind of 
averaging is used in defining the macro quantities, such as concentration, velocity, velocity 
divergence, strain-rate, etc. Second, how to obtain the required constitutive relations for 
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averaged terms like stress, pressure, phase interaction, etc. Neither of these two issues is 
settled at the moment. 

We will refer to the mechanics that governs particle’s motion as the “micromechamcs . 
Such micromechanics includes particle-particle and particle-fluid interactions. In order to 
explicitly formulate constitutive relations, knowledge of micromechamcs is essential. 

It is understood that the formulation of constitutive relation is extremely compli- 
cated, because the micromechanics of fluid-solid interaction is itself not well understood. 
Another reason for this difficulty is less obvious but more significant. That is, up to the 
moment, there has not been a set of governing equations in which all terms are interpreted 
with a consistent averaging method, without such governing equations, and a consistent 
bookkeeping to account for the micromechanics. It is impossible to correctly formulate the 
constitutive relations. The authors believe that this is the source of the recent argument 
about the “solid phase pressure” and related phenomena. A survey of the recent literature 
shows that different interpretations have been given to the solid phase pressure. It has 
been suggested to be equal to the (i) fluid phase pressure (Drew 1976), (u) averaged fluid 
pressure around the surface of a particle (Givler 1987), or (hi) a more sophisticated version 
of (ii) with additional consideration of Brownian forces and bulk viscosity (McTigue et al. 
1986). All of which are intended to apply to an arbitrary flow of a fluid-solid mixture. 
Similarly, the phase interaction (m) in Eqs. (5) and (6) has also been modeled in many 
different ways. Essentially, in all the more recent works this term has been modeled as 

(m)=nh + pV(c) ( 7 ) 

where n is the averaged number of particles per unit volume of the mixture, h is the 
averaged hydrodynamic force per particle. The term p is a source of confusion again. It 
has been equated to the (i) averaged fluid pressure over the fluid-solid interface (Drew 
1983), (ii) fluid phase pressure (McTigue et al. 1986) or (hi) hydrostatic fluid pressure 
(Ahmadi 1987). Again, all the above are suggested for general flow condition. 
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Although the existing two-phase flow models appear to be inconsistent, development 

of the above models has provided a great deal of insights. These insights are essential to 
the work presented here. 

In the following, we will give a step by step derivation of the two-phase flow governing 
equations based on the mass-weighted average defined as 



where c* equals to c for the solid phase and 1 - c for the fluid phase. This averaging 
method has been exploited in the theory of compressible turbulence (Van Driest 1951), 
and applied to granular flow (Ahmadi and Shahinpoor 1983). 

The mass-weighted average of any given quantity is the quantity averaged within that 
phase only. For instance, mass-weighted solid velocity at a point is the average velocity of 
all observed particles passing that point. This average is the easiest one to measure in most 
real flows. Moreover, through using this average, a direct bridge between micromechanics 
and the constitutive relations may be established. Applying this average to the solid phase 
stress, terms such as the solid phase pressure will have a clear meaning. Thus a unique 
definition for these quantities in terms of micromechanics is possible. In addition, a number 
of interesting phenomena that seem to defy the well accepted property of fluids are observed 
when a fluid-solid mixture is viewed as a composition of two separate continuums. 

In the present study, we concentrate on the derivation of governing equations, in- 
terpretation of the averaged terms and the development of constitutive relations from 
micromechanics. Derivation of the actual constitutive relations for various flow conditions 
is beyond our scope and present ability. A few exceptionally simple cases will however be 
studied. In order to complete the mathematical modeling, boundary conditions must be 
derived. Due to the existence of the discrete solid phase, derivation of boundary conditions 
is equally difficult as constitutive relations. This is also left to future development. 
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II. GOVERNING EQUATIONS 


Consider an arbitrary control volume V as shown in Fig. 1. Its surface 5 is the 
control surface- For visual reason this control volume is drawn large compared with the 
solid particle’s size. The derivation is not restricted to this size. The particles me not 
necessarily spherical or uniform either. However, in order to simplify the notations, we 
will discuss the case of uniform spherical particles. The essence of the analysis is captured in 
this simplified case. There has been previous work deriving conservation equations for this 
type of control volume for a fluid-solid mixture (Soo 1981). Nevertheless, such equations 
have not been derived in terms of the mass-weighted average nor the interpretation in 
terms of micromechanics given. 



Fig. 1 A control volume V with control surface S = S, + 5/ and 
interned interface S». 

In general, the flow situation is such that the condition in any given control volume 
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from different observations is a random phenomenon. Only an average behavior may be 
described by a deterministic mathematical model. To allow for cases where flow conditions 
may vary over a length scale comparable with the particle size, the ensemble average over 
many realizations of a control volume smaller than the "representative volume” commonly 
used in fluid mechanics may be necessary. 

For one realization, the rate of increase of solid mass in the control volume is 

dt J v f>' cdV + J s PsCu ndS = 0 ( 9 ) 

where n is the unit outward normal. After ensemble averaging, the resulting integrands 
are smooth functions. Applying Green’s theorem to the second term and removing the 
integrals, the mass conservation of the solid phase in this control volume is obtained as 

dip 9 ) 

-^ + v . ( p -)(„ }= 0 ( 10 ) 


In the above, (p.c) is replaced by (p») and (p.cu) is replaced by (p’){u) where { } is the 
mass-weighted average defined in Eq. (S). Similarly, for the fluid phase, 


dtP) 

dt 


+ v -(^ / ){v} = 0 . 


( 11 ) 


The momentum conservation equations require representation of the forcing terms 

which include the surface and body forces. This conservation for the solid and fluid phases 
in one realization are, respectively, 

dr t 

ftJ v PsCudV + J sPa cuu.ndS = J s cT t .ndS + J^p.c&dV + J^mdV (12) 

dt J v Pf ^ 1 ~ C ) V dV + J s pf ( 1 ~ c ) vv hdS = j( 1 - c)Tf • n dS 

+ J PfO--c)%dV- f mdV (13) 
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Again the ensemble average is first applied to smooth the above integrands. The 
Green’s theorem may then be used to change the surface integral to volume integral. After 
removing the integral sign, the above equations become 

-|b/> s cu) + V • (p 4 cuu) = V • (cT.) + (p«cg) + (m) (14) 

at 

— (pf( 1 - c)v) + V • (pf(l - c)vv) = V • ((1 - c)T f ) + {pf{ 1 - c)g) - (m) (15) 

dt 

Substituting u by {u} + u" and v by {v} + v" and making use of the mass conservation 
equations, we obtain the following equations with the Reynolds stresses for both phases. 
In these equations mass-weighted velocities and stresses appear. 

(p*){d{u}/dt + {u} • V{u}) = V • (c){T s ] - V • (p’){ u"u"} + (p s ) g + (m) (16) 

(p f )(d{v}/dt + {v} • V{v}) = V • (1 - c){Tf} - V • (/){v"v"} + (p f ) g - (m) (17) 

The Reynolds stress in the solid phase is also called the kinetic stress in the granular flow, 
terminology. 

In a realization over a period of time, particles cut by the control surface may in- 
teract with neighboring particles through collisions. The rate of momentum transfer to 
the interior of the control volume resulted from these collisions is part of the surface force 
T s . Moreover, the hydrodynamic forces acting on particles cut by the control surface also 
produce surface forces that contribute to T 8 . We denote these two stresses as T c and T p 
respectively. 

Concept for modeling the mass-weighted average of T c has been described in the 
granular flow literature (e.g. Bagnold 1954, Jenkins and Savage 1983, an excellent survey 
to appear by Campbell 1990). Although most of the work deals with negligible fluid 
effect, the route to extend to fluid-solid mixture is, though complicated, quite clear (e.g. 
Ackermann and Shen 1982, Shen et al. 1988). On the other hand, the explanation of 
the hydrodynamic stress on the solid phase is not readily available. Intuitively, one would 
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quite comfortably accept that the fluid only acts on each particle in the control volume 
through the drag, added mass, etc. Hence, it should contribute to the body force only. 
Consequently the effect of fluid-solid interaction is the hydrodynamic force per particle 
multiply the number of particles per unit mixture volume. This concept is proved doubtful 
as demonstrated by the existing various models discussed in Eq. (7). In the following we 
will rigorously formulate the fluid effect in a two-phase flow. 

Consider a surface particle P in Fig. 2. Part of this particle, P\ is inside the control 
volume, part of it, P° is outside of it. The hydrodynamic force acting on P produces a 
pair of internal forces, ± t, on the intersection of the particle and the control surface S a . 
The total hydrodynamic force acting on this particle, h, can similarly be decomposed into 
that on the outside of P, d°, and that on the inside of P, d‘. 



Fig. 2 Decomposing the hydrodynamic force on a surface particle P. 

The total hydrodynamic force is h = d 1 + d°. 

In budgeting the total force acting on the solid portion of V , the interface force t 
acting on S, most naturally belongs to the solid phase stress. In fact it may equally well 
be classified as part of the body force since after all it acts on P* which is inside of V. 
As long as the budgeting of all forces is done in a consistent way, the resulting equation 
of motion should not depend on the detail of the bookkeeping. Most of the existing 
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two-phase models appear to fail in using a consistent bookkeeping. That is, force acting 
on the control surface and that acting on the internal interface are not always carefully 
distinguished. We choose to call t the surface force for V and d 1 the body force for V due 
to fluid-solid interaction. 

The internal force t is the difference between d° and the total surface force acting on 
P°. The first may be determined if the hydrodynamic force distribution on a particle’s 
surface and the particle’s location relative to the control surface are both given. 

The second is determined using Newton’s second law. That is, the total surface force 
acting on the partial particle P° equals to its inertia subtract the body force acting on 
it. Ensemble average of the inertia and the body force on P° may also be obtained if 
the particle’s relative position with respect to the control surface is given. For a control 
volume reasonably away from the boundary of the flow field, the position of particles on 
the control surface S 8 may be assumed to uniformly distribute inside or outside the control 
volume. Using these arguments it has been shown in Hwang and Shen (1989a) that the 
ensemble average of the sum of t in a unit area produces the following stress, 

{T p } =yXL <e) ■ &rdA - X (v ■ < s » r<iF )' < is > 

where A 0 and V Q are the surface area and volume of particle P respectively, £ is the local 
fluid stress on the particle’s surface, and r is the position vector. The very same equation 
has been obtained by Batchelor (1970) with a quite different derivation. This term has 
been named the “particle-presence stress” in Hwang and Shen (1989a), and it contains the 
“interaction stress” in McTigue et al. (1986). The second term inside the parenthesis in 
Eq. (18) is related to the particle’s rotational inertia (Hwang and Shen 1989a). 

Next we discuss the interaction term m in the solid momentum equation. The most 
natural way consistent with the above bookkeeping is to define it as the total force acting 
on the interface of the two phases inside the control volume V. This interface is denoted 
by Si in Fig. 1, which consists of the surfaces of whole particles if they are entirely in V, 
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otherwise only the part of particle’s surface that is inside V belongs to S',. Total force 
acting on Si is the sum of all h, the hydrodynamic force on a whole particle, and all d‘, 
the portion of the hydrodynamic force on P inside of V. As shown in Hwang and Shen 
(1988), the ensemble average of the total interaction force m is 

( m ) = ^{ h } + <c>V-{ Tf }-V-c{ TP } ( 19 ) 

In the above, {T f } = {Tf} - p/{v"v"} is used to indicate the total fluid stress. Substi- 
tuting Eq. (19) and T s = T c + T p into Eq. (16) , the solid momentum equation becomes 

0 >')(a{u}/a<+{u).V{u}) = V({c){T c })-V( / ) I ){u"u"} + (/)g+M{h} + (c)V{T f } 

r o 

( 20 ) 

This equation is identical to that derived in Prosperetti and Jones (1984). In that work, 
a couple of assumptions were made to arrive at Eq. (20). It is shown here that those 
assumptions are not necessary. Similarly, the fluid momentum equation may be derived as 

</>(dM/^ + {v}.V{v}) = (l — 2 (c))V.{T f } — {T f }.V(c) + (/)g— -^{h}+V-(<c){TP}) 

V o 

(21) 

As it shows in above equations, the effect of particle-presence stress {T p } appears only in 
the fluid momentum equation. 


Derivation of the energy equations follows the same spirit as the above. Hwang and 
Shen (1989c) gave the following equations in indicial form for the turbulence energy in the 
solid and fluid phases respectively, 

fc # >( J + = < c X TC > : v{u} + < c H TP ) ; 

- (p'){u"u"} : V{u} + V • «c){T c ■ u"}) + V • ((c){T p • u"}) 

r 


- V • ({/){— 5— • «"}) + (m • u") - (c){ 7 *} 


(22) 


d , • v" 


2 ^ + ^ 2 V X = “(! - c){pf}V ■ (v) - </){v"v"} : V{v} 


y dt 1 2 

V •((/){( 




J" 




2 p f 

— p(l — c){Vv” : Vv"} — (m • v") 


(23) 
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where 7 * is the sink of solid turbulence energy from particle-particle interaction (e.g. 
collisional dissipation), pf is the fluid pressure, p f" = pf — {p*f}, and p is the dynamic 
viscosity of the fluid. 

The interaction force m does unequal amount of work to the turbulence energy in the 
fluid and solid phases. This difference, as given in Hwang and Shen (1989c), comes from 
the slip velocity between the two phases. 

The combination of mass, momentum and energy conservation equations given by the 
six equations: Eqs. (10), ( 11 ), (20), (21), (22) and (23) governs the six unknowns c, u, v, 
{u" - u"/2}, {v" • v" / 2} and {p^}. To form a closed system, a large number of constitutive 
relations must be obtained. Most of these constitutive relations require knowledge beyond 
the current understanding of fluid mechanics around a particle. One example is the term 
{m • u"} appeared in Eq. ( 22 ). As shown in Eq. (19), m contains the hydrodynamic force 
h acting on a particle. If one considers the drag force part of h, which is in general in 
terms of u" — v", the term {m • u"} will produce {v" • u"}. Namely, the correlation of 
the fluid and solid turbulence must be formulated as part of the constitutive relations. A 
number of analyses have indicated that this correlation depends greatly on the particle’s 
size and density (Chao 1964, Xie 1987, Abou-Arab and Roco 1988). Quantitative study 
of such correlation is far from complete. 

Despite the complexity and the lack of necessary knowledge in fluid mechanics, a 
number of results may be obtained for the constitutive relations appeared in the governing 
equations. In the next section, we will derive a term associated with the fluid-solid in- 
teractions. Namely, the particle-presence stress T p including a detailed discussion on the 
solid phase pressure. A number of interesting observations will be given at the end of the 
derivation. 
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III. CONSTITUTIVE MODELING 


To simplify the notation, we remove all the average symbols ({ ) and { }) and adopt 
indicial notation from this section on. 

In this section, we will discuss the particle- presence stress Tf-. Detail of this term is 
given in Eq. (18). The term in Eq. (18) related to the particle’s rotational inertia cancels 
with the rotational contribution in the particle’s Reynolds stress in Eq. (16) (Babic 1989). 
Hence the net particle-presence stress is only the first term on the right side of Eq. (18). 
This term may be quantified if the distribution of the hydrodynamic stress E is known. 
Before quantifying this stress for special cases, we will provide an interpretation of the 
solid phase pressure first. 

As defined in Continuum Mechanics, pressure means the negative average of the nor- 
mal stresses. Therefore, in addition to the contributions from the collisional and Reynolds 
stress, the particle-presence stress will also produce the following solid phase pressure, 


r2 jt r-n 

I? = -7TT / / Sijtnjtn, sm<f>d<f>d0 

u Vo Jo Jo 

l r 2ir r n 2 

— J J (~P f 6ik - -peuSik + 2^ejfc)n*nj s\n<j>d4>d8 

2 y27T 1 7T 

= 7 -I / p f sin <f> d<f> dd 

4?r Jo Jo 


(24) 


where <f> and 6 are the polar and azimuthal angles respectively. In the above, a Newtonian 
fluid is assumed such that 


£«* = ~P f f>tk - -peuSik + 2 pe ik (25) 

where is the component of local strain-rate. Moreover, for Eq. (24) to be valid, the 
fluid flow near the particle must be incompressible and the strain-rate must possess certain 
symmetry property. Thus the viscous contribution vanishes from Eq. (24). Under these 
special conditions, it is seen that the solid phase pressure is the numerical average of the 
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fluid pressure around the particles. This very result has been postulated by Givler (1987) 
and derived by Prosperetti and Jones (1984) using a couple of assumptions. 

We now derive the particle-presence stress in two extreme cases for a very dilute 
mixture. The two extreme cases are: Stokes regime and inviscid regime. In both cases we 
assume the fluid is Newtonian and the flow is incompressible. 

In case of a shear flow in the Stokes regime, the local strain-rate around a spherical 
particle in an infinite flow field is (Batchelor 1967) 

, N ^ R 5 \ , (rirk&ji + rjr k Su 2 r k r t 5 R 3 5 R* \ 

e 0 (r.) = £„(l - 7 r)+ £*,( ^ 3— ( - *3- + — ) 

1 \/25 . IP 35H 5 \ 


+ Eki 


rki~i 




'A 2 r 3 2r 5 ) ^ 

where r, is the ith component of a position vector r, r =| r |, E t] is the component the 
undisturbed strain-rate and R is the particle’s radius. 

The fluid pressure around a spherical particle in a uniform incoming flow Uoo is 
(Chester and Breach 1969) 

3,. 3 


p = P r + [— 1(1 + |fl,)cos^ + || fl, cos 2 ^ + 0(Rl IogB«)], 


(27) 


where pf is the undisturbed fluid pressure; Uoo is equivalent to the relative velocity of the 
particle to the fluid flow, or U 0 0 = v r where v r =| v — u |; R e is the particle Reynolds 
number defined as 2pfRv r /p\ and <f> is the polar angle of a point on the surface of the 
particle measured from the direction of Uoo- 

Substituting Eqs. (25) and (26) into Eq. (18), the viscous contribution of the particle- 
presence stress is obtained as 

Tff = |(2;i£,j) (28) 

As discussed in Batchelor (1970), this stress together with the viscous stress in the fluid 
phase reproduces Einstein’s formula for the effective viscosity (1956). The solid phase 
pressure is shown below after substituting Eq. (27) into Eq. (24), 

9 , 


f =P f + M PfV ‘ r ' 


(29) 
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In the above, we used p p to indicate the fact that only particle-presence stress T p is 
included in the total solid stress. Neither collisions] nor Reynolds pressure is considered. 

The above is not the whole story about the particle-presence stress. If one by-passes 
the solid phase pressure and investigate the normal stresses directly, the following is found. 


rrPP _ rppp 
x 22 — 1 33 



27 

160 


Pf v r 


(30) 


The total particle-presence stress is the sum of Eqs. (28) and (30). It is apparent that 
Eq. (29) can be reproduced by taking the negative average of Tff, and T-f-f . 

Eq. (30) shows that on top of the viscous effect the isotropic fluid pressure induces 
anisotropic norma] stress in the solid phase stress. This anisotropy of the normal stress is 
the product of the distinctly different dynamics of the two phases involved. If both phases 
move with exactly same velocity, this phenomenon will disappear. 

In case of an inviscid flow, only the fluid pressure contributes to the particle-presence 
stress. This pressure around the spherical particle is (Lamb 1932), 

1 19 

P = P°o + -pfV 2 r (-- + -cos2<l>). (31) 

Substituting Eq. (31) into Eq. (25) and neglecting the viscous part, Eq. (18) becomes 

JR— 

rpp rpp 

1 22 ” J 33 

= ~P f + \pf v l (32) 

and the solid phase pressure from the hydrodynamic effect alone is the negative average 
of the above, 

f = p f - ^p f v 2 (33) 
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The anisotropy of the normal stresses is again observed. Theoretical speaking, without 
considering Reynolds stress, one may observe the shear force in a two-phase mixture even 
though the fluid is inviscid. 

Comparing the coefficient of the solid phase pressure in the above two extreme cases, 
one would expect that as the particle Reynolds number increases from nearly zero to 
nearly infinity, the coefficient in front of the solid phase pressure should vary gradually 
from to — Explicit determination of this coefficient depends on the knowledge of the 
hydrodynamic stress distribution around the particle. Unfortunately this information is 
not available for the entire range of the particle Reynolds number. However, for particle 
Reynolds number greater than 3000, experimental data shows that the viscous contribution 
is negligible. In this case, sufficient information exists to empirically determine the pressure 
distribution, and accordingly the solid phase pressure (Hwang and Shen 1989b). 

IV. SOME INTRIGUING POINTS 

In deriving the governing equations, a few points have struck the authors as being 
quite non-trivial. Some of those may have significant implications that is unclear at the 
moment. 

First, the governing equations are for the mass-weighted quantities. All constitutive 
relations must eventually be written in these quantities to produce a closed system. The 
kinematic quantities appeared in these governing equations are the mass-weighted veloci- 
ties. On the other hand, in constitutive relations, we look for mathematical description of 
the stresses. The average stresses, at least for the viscous component of the fluid stress, 
is however a function of the average local fluid strain-rate instead of the gradient of the 
mass- weighted velocity. Namely, this stress is in terms of the mass- weighted strain-rate in 
stead of the gradient of the mass- weighted velocity. For a fluid-solid mixture with rigid 
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particles, this relation has been derived as (Hwang and Shen 1989c) 

,dvj d{v 3 } 5(1 - c) {u,} Id, . (c) 

W ~ “^T + (i-c) + TT^c ) ^-( <c){u >)) - ( 34 ) 

where is the spin tensor (or particle rotation) of the solid phase. In addition to 

the above two possible definitions for the phase strain-rate, a third strain-rate is also 
important. This is called the bulk (or mixture) strain-rate defined by 


Eij = (1 - c){t f ij) + (c){e'y } 
= (1 -c){e f ij) 


(35) 


where ij and e , ^ are the local strain-rate for the fluid and solid phase respectively. In 
case of rigid particles, e* ; =0. In representing the effective viscous stress in a fluid-solid 
mixture, the above bulk strain-rate must be used. 

In a pure fluid flow, such distinction shown in Eq. (34) or Eq. (35) does not exist. In 
a two-phase flow, this point becomes essential in many places. In fact, if one interprets 
the bulk viscous stress using d{vj}/dx t instead of {dv } /dx t }, the coefficient of effective 
viscosity in Eq. (28), i.e. ^ , will become f or 2 depends on different interpretations of the 
mixture strain-rate. 

Second, in Givler (1987), a phenomenon called anti-diffusion in the solid phase was 
mentioned. This term is believed to come from an inconsistent bookkeeping system of the 
surface and interaction forces in the control volume V. The equation given there for the 
solid phase was 


P \au,/9t + = p'm + (k - - & - ?') 


(36) 


When the potential flow is considered, the dc/dn term has a positive coefficient. 
This presents a force that moves the solid phase from a low concentration towards high 
concentration, which is a puzzling phenomenon. However, in the present derivation, it 
is shown that in the solid momentum equation such dc/dxi term does not exist. The 


50 



associated diffusive or anti-diffusive forces are thus absent. But, investigating the fluid 
momentum equation, we do observe such dc/dxi term. Through the mass conservation, 
any diffusion in the fluid phase results an effective anti- diffusion in the solid phase and vice 
versa. At this point, we have observed a totally opposite trend between Givler (1987) and 
the present model. Namely, Givler’s model would predict a diffusion phenomenon for solid 
particles in the Stokes regime and anti-diffusion phenomenon in inviscid flow regime. The 
present model however predicts exactly the opposite. Further verification with detailed 
experimental data is desirable. 

On top of the above, the anisotropy of the solid phase normal stress produced by fluid 
pressure is a new observation. This may produce interesting thermodynamic properties 
that axe peculiar to a two-phase flow only. 

It is natural to ask whether such a great care in deriving governing equations shown 
here is of importance. In order to see this, two models have been applied to a vertical pipe 
flow of a fluid-solid mixture (Hwang and Shen 1988). The two models axe identical except 
the phase interaction term m,. For case (A), m; is modeled as shown in Eq. (7) with p = pP 
and case (B) is the present model given by Eq. (19). The resulting non-dimensional slip 
velocity u* — v* verses the non-dimensional fluid pressure gradient (^f )* is reproduced in 
Fig. 3 for three different density ratios p* = &-. Qualitatively different results axe obtained 
from these two models. From the behavior of the neutrally buoyant particles, it is believed 
that the present model is more reasonable. 

V. CONCLUSION 

For a two-phase flow, the interpretation of terms in the governing equations is the first 
step towards deriving constitutive relations using micromechanics. Therefore it has been 
surprising to us that the work presented here has not been available in the vast amount of 
two-phase flow literature. On the other hand, this may be explained by the fact that only 
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recently researchers have started to derive the constitutive relations from micromechanics. 
Hence the need for such derivation of governing equations is also very recent. 



Fig. 3 Effects of p,* and ( dp/dx )* on slip velocity from 
model A and model B. 

It is shown in this work that the derivation of governing equations for flows of a fluid- 
solid mixture is not as straightforward as any single phase continuum. This derivation 
is done with a careful and almost tedious method to account for all transfer quantities 
between each phase. The resulting equations in terms of mass-weighted average quantities 
differ from the existing models. 

These governing equations provide a foundation to incorporate micromechanics in 
deriving the constitutive relations. A number of previously well accepted facts about fluid- 
solid flows emerge naturally from this analysis. In addition, interesting phenomena such 
as: non-equality of average strain-rate and gradient of the average velocity, anti-diffusion, 
and anisotropic normal stress have been observed in these equations. Interpretation and 


52 



implications of these phenomena is of great interest. 
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STRESS IN DILUTE SUSPENSIONS 

Stephen L. Passman* 

Abstract. Generally, two types of theory are used to describe the field equations for suspensions. The 
so-called “postulated” equations are based on the kinetic theory of mixtures, which logically ought to give 
reasonable equations for solutions. The basis for the use of such theory for suspensions is tenuous, though 
it at least gives a logical path for mathematical arguments. It has the disadvantage that it leads to a 
system of equations which is underdetermined, in a sense that can be made precise. On the other hand, 
the so-called “averaging” theory starts with a determined system, but the very process of averaging renders 
the resulting system underdetermined. I suggest yet a third type of theory. Here, the kinetic theory of 
gases is used to motivate continuum equations for the suspended particles. This entails an interpretation 
of the stress in the particles that is different from the usual one. Classical theory is used to describe the 
motion of the suspending medium. The result is a determined system for a dilute suspension. Extension 
of the theory to more concentrated systems is discussed. 


1. Introduction. In theories of multiphase flows, it is natural to postulate or to 
derive equations of balance similar to those occurring in the theory of dilute mixtures 
of gases [1,2]. The usual process of doing so, along with reasonable assumptions for the 
constitutive properties of the materials composing the flow, always leads to a system with 
more unknowns than equations. Though there is no definitive reason that this is a bad 
situation, intuition abetted with proved theorems for special types of systems indicates 
that the normal desirable situation is the same number of equations as unknowns. The 
resulting quandary for multiphase flows is known as the “closure problem”, and methods 
for “solving” or “closing” it, that is, finding “sufficient” additional equations, has been 
the focus of considerable research in multiphase flows. Here, we try to shed some light on 
such problems. Essential to doing so is stating the problems unequivocally. In order to do 
that, we choose a special but interesting physical situation, then give typical equations of 
balance and constitutive equations for that physical situation, according to a continuum 
theory and an averaging theory. The closure problem occurs in both types of theories, 
though its form is different. However, it is possible to formulate theories in which the 
closure problem does not occur, and therefore need not be solved. A physical basis for 
such a system is presented, and a putative set of field equations is suggested. 

2. Determined, Underdetermined, and Overdetermined Systems of Equa- 
tions. Assume we have a system of equations of the form 

ft ( y j 5 Dkyj) = o, 

with i = 1 , j — and k = 1 The fi are n functions of the m 

variables yj and their derivatives up to order p . The system is called determined if n = m, 

^Pittsburgh Energy Technology Center, Pittsburgh PA 15236, on temporary assignment from Sandia 
National Laboratories, Albuquerque NM 87185. 
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overdetermined if n > m, and underdetermined if n < m. By our definition, all systems 
of equations are of one of these three types. Ideally, of course, it would be convenient if 
determined systems always had solutions and they were unique, if overdetermined systems 
never had solutions, and if underdetermined systems always had families of solutions. That 
this is not the case can be shown by examples. To begin, consider the underdetermined 
system 

(1) x\ + x\ = 0. 

Naturally, specifying a system of equations is meaningless without specifying their domain, 
but since this paper is informal, I follow the convention of doing so tacitly, that is, all 
functions are mappings of all real variables for which they can be defined reasonably into 
a range defined by the function. Here, the underdetermined algebraic system (1) has the 
single unique solution 

X\ =0, X2 — 0. 

Now consider the determined system 

^1+^2 = 1, 

2x^ T 2x2 = 2. 

This system does not have a unique solution, rather it has an infinite one-parameter family 
of solutions. Finally, consider the overdetermined system 

X\ + x 2 = 1, 

2xj -f- 2 x 2 — 2, 

2xi T 4 x 2 = 4. 

This system has a unique solution. All of the examples cited involve algebraic equations, 
not differential equations, but of course examples of the same type can be constructed with 
differential systems. 

It is easy to object to the arguments above because the systems cited are “special” or 
“pathological”. Indeed, I agree with that type of objection, and in a sense that is just the 
point of this discussion, for to make such arguments, one must use very special properties 
of the systems. Furthermore, for each system, it would have been possible to rearrange the 
system using simple manipulations, obtaining very complicated new systems with exactly 
the same properties . Proofs of existence or non-existence, uniqueness or non-uniqueness, 
then would be much more elaborate exercises, perhaps depending on sophisticated math- 
ematics or luck for ultimate outcome. 

A different line of argument is possible. Most often, the equations governing multi- 
phase flow are systems of partial differential equations, so complicated that they are not 


( 3 ) 


( 2 ) 


-i» J 
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easily amenable to existence or uniqueness theorems. Such equations may be relevant to 
important problems in technology, so important that they must be solved, in no matter 
in how vague a sense, immediately. More often than not, that means the use of large 
computer codes. Despite the above examples, we find that an intrinsic part of building 
such codes is the desire of the numerical analyst for determined systems, that is, systems 
determined in exactly the sense defined here . The idea that the equations obtained by 
specialists in multiphase flow are “independent” is supported in some vague sense by the 
fact that they have different physical meanings. For example, some are balance equations 
for the constituents, some are constitutive equations, and some are constraints. 

3. Continuum Theories of Multiphase Flows, Here, we consider a standard type 
of theory for multiphase flows, as derived from continuum considerations. 1 Intrinsic to such 
considerations is the assertion that each constituent fills all of a region of space. This is the 
basic assumption of theories of interpenetrating continua or “solutions” [4]. The theory 
then is made to model a multiphase medium by the inclusion of volume fractions <f) a as 
basic variables. A typical set of field equations for such a continuum having n constituents 
is 


( 4 ) 


£"=l^a =1, 

V 

<f>a + 4>a div V a = 0, 

PaV a = Pa^a + m a + divT a , 

£” =1 m a = 0, 

{T a ,m a } — g(v\>,(})b, and their derivatives). 


Here for simplicity we consider only pure mechanics, and the symbols have the obvious 
meanings. The first equation expresses the fact that the material is saturated; the second 
and third are balances of mass and momentum for each of the constituents. Each con- 
stituent is assumed to be incompressible, and the p a are the reactions to those constraints. 
The fourth equation is conservation of mass for the mixture, and the last equations express 
constitutive properties of the constituents, in particular the dependence of the stress on the 
deformation rate and other properties for each constituent, and appropriate expressions 
for interactions of the two materials. We note that these last expressions can be somewhat 
problematical, and in fact debate about their forms has generated a considerable litera- 
ture. They are not discussed here. Rather, we assume they are known for applications of 
particular interest. 2 Our intuitive feeling from considering the physics of this situation is 
that such a system of equations is “complete”, but in fact that is not the case. Here and 
henceforth, for the purpose of counting equations, we assume the multiphase flow consists 
of two constituents. Though such is not always the case, for the purpose of our argu- 
ments here, that case is general. The result is a system of 9 equations in the 10 unknowns 

1 Discussion of the basis of such theories, as well as references to the standard works, are given in [3]. 

2 See [5] and [6] for a discussion of these equations. 
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{pat </»fl,v a }, that is, the system is underdetermined. The usual physical motivation for this 
apparent quandary is plausible: Though Equations (4) express the exchange of momentum 
between constituents, that is not the only way the constituents interact, for in addition , 
there should be a force balance between the constituents. The most primitive visualization 
of this is sufficient for arguments here. That is, the solid phase is considered to consist of 
spherical particles of one size, surrounded by the fluid phase. Then a radial force balance 
on a single particle gives 

(5) P*=Pf ■ 

This obvious and elegant closure argument gives a system of equations which is notori- 
ously ill-behaved [7,8], so much so that it must be rejected. More sophisticated arguments 
are possible, and they sometimes appear to suffice to render the system of equations thus 
obtained to be at least well-behaved enough to be handled by standard computational 
techniques. Usually, the arguments adduced are generalizations of those leading to Equa- 
tion (5) in that they consider a particle in a flow field of a known type at infinity, then 
use techniques of hydrodynamics to solve or partially solve for the flow field around the 
spherical particle. Surface tension may be considered also. Of course the resulting pressure 
on the particle is a function of position on its surface relative to the flow field at infinity, 
so some sort of integration is required. The result is 

( 6 ) P» =Pf + /(n,a), 

where II denotes properties of the flow field, and a denotes properties sufficient to charac- 
terize surface tension. I note that since / depends upon the flow field at infinity, adducing 
it as a constitutive relation valid for all flows has the potential for leading to inaccurate 
results. 3 

In addition to the mathematical argument against using Equation (5) as a closure 
relation, there is a physical argument against it, which also is inherited by Equation (6). 
In doing the arguments leading to these equations, the assumption is that p s and pf are 
pressures “in” the respective materials, and that it is appropriate to write an expression 
for one in terms of the other. The need for the closure relation comes from arguments 
about the system (4). In this system, the pressures are derived as reactions to constraints. 
Therefore [10], they are dependent variables of the system of equations, totally independent 
of one another. In treating the complete system of equations and boundary conditions. 4 
the quantities p s and pj thus cannot be related a priori . Another way to see this is that in 
fact Equations (4) are field equations for the whole continuum , while the closure relation 

3 A similar difficulty arises in rheology, where it is commonly known [9] that no constitutive equation 
giving an accurate representation of the physics of shearing flows also represents stretching flows adequately. 

4 Boundary conditions in themselves constitute a difficult problem for multiphase flows. They are not 
discussed in this paper. 


60 



(6) is not. Rather, it is derived from a “micromechanical” argument, then scaled up in a 
way the nature of which never is made clear, so that the symbols p s and pj in Equations 

(6) are assumed to have the same meaning as the same symbols in the closure relation, 
without proof or explanation. 

4. Averaging Theory. The basic ideas behind averaging theories are diametrically 
opposite from that of the continuum theories, though the objective — finding differential 
equations for fields, valid throughout a body — is exactly the same. It is perhaps fortuitous, 
or perhaps a sign that the equations actually represent some sort of physical “truth”, that 
the forms of the equations resulting from the two approaches are so similar. For averaging 
theories, the region of space occupied by the material is thought of as being occupied by 
two different types of body, the suspended particles and the suspending medium. Each of 
the types of body is considered to be distinct in the sense that the join [11] of the bodies 
constitutes all of the space occupied by the composite body, while the meet is empty. Then 
each of the types of body is an ordinary continuum, and satisfies exactly the balance and 
constitutive laws expected of an ordinary continuum, that is, 

divv a = 0, 

(7) p a V a = pab a + divT a , 

{T a } = g(vb, and their derivatives). 

Here, of course, the bodies still are capable of momentum interaction, but unlike the previ- 
ous situation, the micromechanical model for momentum interaction has a clear meaning. 
This is eight equations in eight unknowns, and thus is a determined system. Moreover, 
conditions for the difference of pressure such as Equation (6) now have a correct theoretical 
status, for now they are not field equations, rather, they are boundary conditions. Thus a 
determined system is obtained, and it is mathematically correct and physically plausible. 
The difficulty, of course, is that to formulate a boundary-value problem, a reasonable set 
of boundary and initial values for every particle in the system is needed. Such information 
normally is not available for any physical problem. Even if it were, finding a solution, with 
or without a computer as an intermediary, would be a nearly hopeless task. Moreover, 
even if such a solution were found, most of the information it contained would be of little 
use, because it would be too detailed. A plausible way to digest such data would be to 
average it in some sense. The usual approach in averaging theory is, not to go to the 
considerable trouble of averaging the solutions to (7), but rather to average (7) and then 
solve the averaged equations. Such an approach is highly appealing intellectually, but is 
fraught with mathematical difficulty. This paper is not the place to discuss such difficulties 
in detail. One, of course, is that the term “average” which has been used in a very vague 

5 Of course, the same argument can be made for the expressions in (4) for m a . There, however, the status 
of the equations is clear, because the appropriate micromechanical arguments can be used as motivation 
for the continuum theory, which then gives exact relations having the same status as field equations. 
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way here, must be given a precise meaning. It is fortuitous that, for most of the averaging 
methods tried so far, the resulting equations have almost the form of the Equations (4) 
derived from the continuum theory. Unfortunately, for every method of averaging I have 
seen, though a determined system is averaged, the result of the averaging process is an un- 
derdetermined system, that is, the averaging process makes the closure problem reappear. 
Most readers will be familiar with why this happens without going through the details, 
for the averaging process always is similar to that used in turbulence theory and some of 
the extra terms are of the same form as Reynolds stresses. In a broad sense, then, though 
the continuum theory and the averaging theory start from different places and proceed by 
different methods, they end in approximately the same place: underdetermined systems of 
approximately the same form. 


5. Sketch of a Theory for Dilute Suspensions. Previously in this paper, much 
has been made of the fact that most of the equations in the continuum theory have been 
“postulated”. It is possible to interpret that terminology as meaning that they have been 
made up with no mathematical or physical basis. In fact, that is far from true. The 
kinetic theory of dilute monatomic gases for identical gas molecules is well-known and 
is commonly taught in courses for graduate students in science and engineering [1,2,12]. 
Much less well known is the fact that there is a similar theory for gases where there are 
a finite number of different types of molecules; in other words, a solution of several gases 
[13,14]. The resulting balance equations are exactly identical to those for the postulated 
theory of mixtures. 

Here, I use the motivation of the kinetic theory of gases to support a mixture the- 
ory in an entirely different way. Most important for the discussion here is the fact that 
there is an exact definition of the pressure, and it is not the pressure “in” the particles, 
rather it is a momentum flux — an entirely different concept. 6 Moreover, it is possible to 
force agreement of the theory with that of a viscous compressible gas, with the viscosity 
determined in terms of molecular parameters. I consider a dilute solution of particles in 
an inviscid fluid. Consider only the particles. They axe an agglomeration of molecules, 
exactly like those in the theory of a monatomic gas, except that the scale of the molecules 
is somewhat larger than in a gas. Thus, precisely the same arguments can be used to 
motivate a continuum theory for the particle phase of the multiphase flow as is used for 
a gas. All of the expressions are the same, and e.g ., one can accept the viscosity of the 
particle phase as a phenomenological coefficient, or one can consider it to be determined 
from molecular quantities, according to one’s taste. In either case, unlike in the theories 
discussed in the previous two sections of this paper, it does have meaning. The equations 


6 In another paper in this volume, O. Walton uses the same definition in his computer molecular dy- 
namics simulations. 
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for the particle phase then are 


Ps + p 3 d\v\ s = 0, 
psV s = Psb s + m + div T a , 

(8) 

Ps = Ps(ps), 

T, = T 3 (sym gradv s ). 

This is a system of 5 equations in 5 unknowns, that is, a determined system. Now let the 
molecules be submerged in an incompressible fluid. Naturally, there will be an interaction 
between the particles and the fluid, and this interaction can be expressed as a constitutive 
equation for m, which can be thought of as a part of the body force b/. Of course, the 
equations for the fluid phase are the expected ones, 

div v/ = 0, 

(9) pf\f = p/bf - m + divT/ 

T/ = T/(sym grad V/). 

again, a determined system. Thus, for a theory of this type, no closure problem exists. 7 

Generally in a theory of this type, one expects to see volume fractions appear intrin- 
sically. Since the ideas here are for a very dilute, saturated suspension, the concept is not 
very important, except, perhaps, in the constitutive equation for m and in formulas for 
the “effective viscosity” [5,12], Of course the idea can be introduced formally by setting 

Pf = 7 


with 7 / a constant, and 


4>f + (f> 3 — i. 


These substitutions introduce the same number of equations as unknowns. 
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7 In another paper in this volume, J. T. Jenkins and D.F. McTigue present a similar set of ideas for a 
very concentrated suspension. The kinetic theory is worked out in considerable detail. D.F. McTigue also 
has presented a similar set of ideas in an unpublished lecture at the meeting of the American Geophysical 
Union held in December, 1988. 
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THEORY OF DROPLET: I 

RENORMALIZED LAWS OF DROPLET 
VAPORIZATION IN NON-DILUTE SPRAYS 
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ABSTRACT 

The vaporization of a droplet, interacting with its neighbors in a non- 
dilute spray environment is examined as well as a vaporization scaling law 
established on the basis of a recently developed theory of renormalized 
droplet. 1 The interacting droplet consists of a centrally located droplet and 
its vapor bubble which is surrounded by a cloud of droplets. The distribution 
of the droplets and the size of cloud are characterized by a pair-distribution 
function. The vaporization of a droplet is retarded by the collective thermal 
quenching, vapor concentration accumulated in outer sphere, and by the limited 
percolative passages for mass, momentum and energy fluxes. The retardation is 
scaled by the local collective interaction parameters; group combustion number 
of renormalized droplet, droplet spacing, renormalization number and the local 
ambient conditions. The numerical results of a selected case study reveal 
that the vaporization correction factor falls from unity monotonically as the 
group combustion number increases, and saturation is likely to occur when the 
group combustion number reaches 35-40 with interdroplet spacing of 7.5 
diameters and the environment temperature of 500 K. The scaling law suggests 
that dense sprays can be classified into: (1) a "Diffusively Dense" cloud 

characterized by uniform thermal quenching in the cloud, (2) a Stratified 
Dense" cloud characterized by a radial stratification in temperature by the 
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differential thermal quenching of the cloud, or (3) a "Sharply Dense" cloud 
marked by fine structure in the quasi-droplet cloud and the corresponding 
variation in the correction factor due to the variation in the topological 
structure of the cloud characterized by pair-distribution function of quasi- 
droplets . 


NOMENCLATURE 


D Mass diffusion coefficient 

g Pair - distribution function 

G Group combustion number 

L Latent Heat of vaporization 

m Vaporization rate 

n Number density 

q Heat of combustion 

r Radial coordinate 

R Radius 

s Droplet separation 

T Temperature 

u Velocity 

W Molecular weight 

a Schvab-Zeldovich Variable 

8 Renormalization number r\ /n 

ts co 

y F L /q 

£ f ~ (w f v f )_1 

0 ^ (° ) 

v Stoichiometric coefficient 

n r/r 

i 

p Density 


66 



Subscript 


c Canonical bubble 

co Canonical bubble of test droplet 

F Fuel 

T Temperature 

l Liquid drop surface 

ts Transition sphere surface 
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1. Introduction 


Increasing theoretical findings 1 ^ and experimental evidence 14-18 

attest the widely held belief that the short-range collective interaction 
11 — 13 

among the neighboring droplets and the long-range interaction ^ 9^ 
with the droplets at distance, on a hydrodynamic scale, have profound 
impact on the state of a droplet, i.e. the states of saturation, vaporization, 
ignition, combustion and extinction, as well as the droplet interfacial 
process rates in non-dilute cloud or spray environments* These collective 
interactions, produced by local hydrodynamic and transport processes in a 
complex topological environment, serve to control the percolative passage for 
dispersing mass, momentum energy fluxes and the effective interfacial area for 
the property exchange processes* These interactions result in the collective 
thermal quenching, the accumulation of vaporizing species and the tendency for 
stagnating the microscale local Stefan flow and mean flow through dynamic 
equilibration between the two phases* 

Review of the current theories of collective interactions including: 
Group Vaporization and Combustion l > 7 » 10 (GVC), Discrete Droplet Model ( DDM) 

* ^ arid Droplet In Bubble (DIB) ^ ® reveals two major theoretical 

deficiencies in the theory of short-range interaction and droplet rate 
processes at this juncture. These are: (1) the lack of the fundamental 

concepts and the mechanisms interlinking small scale discrete droplet 
processes with a large scale quasi-continuum flow of a non-dilute cloud or a 
spray, and (2) the incomplete understanding of complex interfacial processes, 
finite rate reaction, turbulence and transport processes that occur in the 
vicinity of each droplet. The first issue compels the current research of 
non-dilute sprays to proceed with two rival approaches; GVC, that primarily 
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deals with the long-range interaction, thus neglecting small scale resolution, 
and DDM, which focuses attention on the detailed behavior of discrete droplets 
in well-ordered droplet assemblies, thereby requiring a laborious compu- 
ational procedure involving a large number of droplets encountered in 
practical sprays. 

The alternative approach, adopting a continuous spray model supplemented 
by realistic droplet laws, appears the compromised mean of the prediction of 
non-dilute spray. A recent study of Tishkoff ^ on the numerically correlated 
vaporization correction factor derived from the DIB model demonstrates the 
viability of such combined approach to complement modern non-dilute spray 
calculations formulated on either Eulerian— Eulerian or Eulerian Lagrangian 
framework. It must be mentioned, however, that an attempt to deduce an 
improved droplet vaporization law from the results of existing DDM met with 
difficulties due to the geniune lack of a self-consistent criterion of 
"droplet environment" whereupon the gas properties such as the temperature and 
concentration of vapor species, are inserted in the droplet law for the 
determination of the vaporization rates. A universal theory of short-range 
interaction and a set of comprehensive laws of droplet rate processes remain 
as the major unsolved issues in the contemporary theory of non-dilute sprays 
which has been the central theme of research 12, 13 conducted. The objectives 
of this paper are to present the basic concepts, theoretical approaches and 
the results of Renormalized Droplet (RND) theory to establish vaporization 
laws for droplets In a stationary non - dilute cloud environment. 

The paper begins with the descriptions of droplet models and theories, in 
section 2, to clarify basic concepts and definitions of "droplet, adopted in 
modern spray theory. The structure, model and mathematical analysis of RND 
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are described in Section 3. The results of numerically predicted RND 

structure, vaporization correction factor and their dependence on collective 
interactive parameters are presented in Section 4. 

2. Droplet Models: Topological Properties and Modeling 

2* 1 Elemental Model 

Three fundamental topological properties that play basic roles in 
interaction phenomena between the droplet under investigation, termed the 
"test droplet", and its neighbors, i*e. M field droplets," are; (1) 

graininess or discreteness" of a droplet, (2) "covolume" or "evacuated 
volume" of a droplet, and (3) "localization" of a droplet relative to a 
reference point. Basic droplet models which attempt to capture properties of 
the droplet to a desired level of accuracy, are classified into three types 
according to the level of sophistication in the characterization of the 
topological features, as summarized in Table 1. First, the "natural droplet 
model which characterizes a droplet by its true geometrical shape, size and 
location, simulates all the topological properties described above. 
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Table 1 
Droplet Models 


Droplet 

Properties 

Natural 
Drop Model 

Quasi-Droplet Mode 

*1 

Non-Unif ormly 
Smeared Drop 

Uniformly 
Sheared Drop 

Point Source 
Model 

Graininess 

Sharp* 

Diffuse** 

No 

Singular 

Covolume 

Sharp 

No“" 

No 

No 

Localization 

Sharp 

Diffuse^ 

No _ 

Singular 


*Sharp properties caused by phase discontinuity 
**Dif f use properties featured by the absence of phase discontinuity 


The natural droplet model provides a sharp droplet configuration required 
for the predictions of microscale flow structure, and results from a 
combination of convective and Stefan flows, and interfacial process rates of 
the test droplet. The solution is given by solving Dirichelt boundary value 
problems associated with conservation equations and a well defined set of 
physical boundaries formed by the interface of droplets. 

The second type of model, i.e. "quasi-smeared droplet," describes a 
droplet by a medium spreading through the space. The droplet thus coexists 
with the host medium without apparent phase discontinuity. The models in this 
category are further divided into "uniformly" and "non-unif ormly" smeared 
droplets, Fig. 1. In the latter model,, the location and spatial extension of 
droplets are depicted by a joint probability distribution function. The 
"graininess" and "localization" of droplets are therefore partially 
characterized by the joint probability function rather than the phase 
discontinuities as in the case of the natural drop model. 

"The uniformly smeared droplet model" which has been adopted in the 
majority of two-phase flow and spray theories is a special case with a uniform 
distribution of the joint probability with its numerical value equals to unity 
throughout the space. In this model, "graininess" and "localization," are 
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absent. Because of the lack of the basic topological properties, these quasi- 
droplet models are unable to be used for the determination of the droplet 
bound flow field and droplet laws. However, when appropriate droplet laws are 
provided to describe the interfacial processes rates, the flow structure of a 
farfield can be predicted with an acceptable accuracy by a quasi-two phase 
approach that is simpler than the Dirichelt boundary value problems. 

The third model is a "point source model" which uses a mathematical 
singularity to characterize the "location" and a "singular graininess" of a 
droplet in a dilute spray with a large spacing, i.e., the ratio of spacing to 
the droplet size is much greater than unity. Like a "quasi-droplet model", 
the point source model" fails to provide a satisfactory near-field structure 
required for the prediction of the rates of interfacial processes, but the 
model offers a simple mathematical theory for the prediction of far field 
solutions. 

2.1.2 Composite Model 

Although a specific elemental model has been frequently used in the 
analysis of single droplet, or many-droplet problems (e.g. natural droplet 
model in DDM and DIB, uniformly smeared droplet in GVC), the potential 
advantages of the simultaneous use of two or more than two elemental models 
has not been fully exploited. 

A composite modeling technique permits a strategic selection of a desired 
elemental model in a certain selected region and an alternative model in 
another region to enhance the modeling flexibility and the simplicity in 
analysis. The choice of an elemental model is determined by the type of the 
flow field data and accuracy required in each region. For example in RND 
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theory, the test droplet is modeled by a "natural droplet model" to facilitate 
the determination of the near flow field, which depends on the detailed 
geometry of the test droplet, and the flow disturbance created by the 
neighboring field droplets. On the other hand, "field droplets" are modeled 
by "non-uni for mly smeared droplet" which are distributed in the neighborhood 
of the test droplet. The distribution of these quasi— droplets is described by 
a joint probability function to simulate localization and graininess of the 
quasi-droplets. This composite model provides a simple, self-consistent and 
useful analytical method of treating the interfacial phenomena of a test 
droplet interacting with its neighbors. The details of the application of 
composite model in RND is described in the remaining part of this section and 
the next section. 

2*2 Theories of Droplet: Canonical and Renormalized Representation of a 

Droplet Under Short Range Interaction. 

"The Theory of Droplet" concerns itself with interfacial phenomena and 
the process rates of a test droplet in isolation or under the influence of the 
collective interaction of field droplets* A representative analytical 
approach for the former is the single droplet theory, and for the latter case, 
the approaches are DIB and RND models* Since the latter problem concerns 
collective interaction between a test droplet and field droplets, the nature 
of the problems and the theoretical procedures are, in general, similar to 
that of "many droplet problems". However, an important distinction between 
the two theories is that the "droplet theory" aims to establish the rate of an 
interfacial process of a test droplet on the framework of "the minimum-sized 
many-droplet system" that includes the smallest number, N m £ n , of droplets 
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i.e., a test droplet and N in ~l field droplets confined in, "the minimum sized 
region", i.e. R, llin region. The requirements of "the minimum number of 
droplets" in "the minimum region" are imposed in RND theory so that the 
droplet laws, deduced from such a minimum droplet system, can be explicitly 
formulated in terms of the self-consistently defined "local properties" of a 
spray flow field. This basic feature of "local representation," is essential 
in modern spray calculations that use local droplet laws. In contrast to the 
theory of renormalized droplet, the "cluster" model 8 *** treats a droplet 
system with an arbitrary number of droplets in a finite region. Thus, the 
model is appropriate for dealing with a partial or complete domain of a spray 
by an approach that is different from the conventional spray theory. 

2.2.1 Canonical Droplet Theory 4 : CDT 

The DIB 4 which consists of a natural droplet and its Wigner-Seitz bubble, 
has two fundamental properties. Firstly, a high degree of symmetry preserved 
in the DIB's periodical droplet assembly inhibits the transports of mass, 
momentum and energy transports among the neighbor droplets. Thus, DIB 
portrays a thermo-chemically closed adiabatic system with respect to neighbor 
droplets. Secondly, the strength of the interaction between the droplet and 
its bubble is determined by the size of the bubble which has the radius equals 
to half of the droplet spacing. When the droplet separation becomes 
infinitely large, the predicted vaporization rate approaches to that of an 
isolated droplet, as expected. 

Because of these two unique features of the DIB which may be regarded as 
a reference model of a short - range interaction, the theory will be termed 
"canonical droplet theory" (CDT). The theory provides a basis of RND model 
described in the following subsection. 
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2*2*2 Renormalized Droplet Theory 

Since the canonical theory treats the interfacial processes of a droplet 
in a hypothetical, well ordered droplet lattice, the validity of the canonical 
droplet law is questionable when applied to practical sprays with disordered 
droplet distributions. Needless to say, the adiabaticity and closure of the 
test droplet relative to its neighbor also break down in practical sprays. 
Another theoretical deficiency of DIB is the lack of the geometrical 

compatibility of the "environment" with that of droplet environment in spray 
theory. In CDT, the edge of a bubble has been adopted as the environment, 
though such choice is not necessarily the unique alternative, and the gas 
temperature or species concent ration is used for the determination of the 
vaporization rate and droplet transient heating rate. However, in a spray 
calculation, the local average gas phase temperature of multi-phase flow is 

used for the determination of droplet process rates. Clearly, the environment 
of a canonical droplet does not coincide with that of a spray. Such 

incompatibility prohibits the encoding of the droplet law, derived by DIB, 
into a spray calculation unless an explicit link between two environmental 

properties, i.e., DIB vs smeared droplet model adopted in a conventional 
multi-phase continuum spray is provided. 

The RND theory provides a rigorous theoretical procedure that removes the 
major shortcomings of the canonical theory and provides encodable droplet laws 
by the applications of (1) composite droplet modeling technique and (2) the 
minimum sized many-droplet system, as described in the following section. 
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3. Renormalized Droplet (RND) 

3. 1 Structure of Model 

The theory of renormalization^ portrays a test droplet interacting with 
its neighboring field droplets by a composite "dressed or clothed" droplet 
model shown in Fig. 2. Two principal structural elements of RND are (1) a 
Dr op let -In-Bubble (DIB) and (2) a cloud of non-uniformly smeared droplets 
which functions as an external clothing for the DIB. The distribution of the 
quasi-droplet is described by the pair-distribution function representing the 
joint probabiltiy of finding a test droplet and its neighbor at a separation 
s. The pair-distribution function vanishes in the immediate vincinity of the 
bubble, representing the evacuated volume effect of droplets, and then 
increases rather rapidly to a value greater than unity at the radial distance 
comparable to a mean droplet spacing. This first high droplet density region 
populated by the nearest neighbor droplets is termed the "first coordination 
shell." The population peaks of the "second and higher order coordination 
shells" lose their sharpness as they merge with one another and are ultimately 
lost in a continuous environment where the pair distribution function 
approaches to unity. The size of the transition sphere, R ts , defined as the 
radial location where the pair distribution function is 0.99, depends on the 
droplet size, spacing and arrangement. The ratio of the size of the 
transition sphere, R , to c ^ e characteristic hydrodynamic scale, L, is 
typically much smaller than unity, and thus the sphere degenerates into a 
point in the limit when R /L vanishes. Accordingly, the average properties 
of the gas over the transition sphere approach to the local properties of the 
continuum flow as the limit R /L goes to zero. This "correspondence 
hypothesis" and the interlinking transition sphere constitute the two key 
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factors for the determination of the droplet process rate by the local gas 
properties of the continuum flow field. In the present analysis, RND is 
assumed to be spherically symmetric, and non-reacting. Transport properties 
are constant and the Lewis number is unity. 

The criterion of a quasi-steady state is a much more complex issue than 
that of a single droplet because of the multi-scale and multi-time phenomena 
linked with mass and energy transport in typical RND. In general, the quasi- 
steady state assumption is valid when (1) the largest characteristic diffusion 
and conduction times associated with the transition sphere is much smaller 
than the life time of any droplet in the cloud, (2) no droplet in the 
transition sphere is in the state of transient heating, and (3) no gas phase 
region is experiencing an initial or a terminal transient process. The 
effects of transient processes, and the validity or the limitation of the 
quasi-steady theory in a practical spray are discussed in Section 5. Detailed 
time-dependent analysis of transient processes in RND, which will be presented 
as Part II in the future, reveal that RND is found in some dense sprays to 
exhibit only a brief or finite period of quasi-equilibrium state. 

Additionally, RND exhibits dynamic saturation in time scales comparable to the 
characteristic diffusion time in the canonical bubble when the local group 
combustion number exceeds a critical value. 


3.2 Mathematical Analysis 

Non-dimensional equations governing RND are formulated separately for the 
inner DIB region and an external quasi -droplet cloud, respectively, as 


follows. 
1 d 


2 d h. 
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where n = r/r^(0), p and v are non-dimensionalized by the gas density, and 

3 

velocity on the droplet surface, respectively, Gs = 4*imr , g is a pair 

v.o 

distribution function, p is the vaporization shape factor defined by (12), 
t are the fluxes of various properties; mass for i = m, fuel vapor for i=F 
and thermal energy for i=T, and A refers to the properties pertaining to an 
outer region. The definitions of properties ou and constants are 

summarized in Table 2. 
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Schvab-Zeldovichs Variables and Constants 
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The system of Eqs. ( 1 ) — ( 4 ) is integrated by a repeated quadrature. The 
constants of integration are determined from the conditions of the 
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impermeability of the gaseous species on the droplet surface and the balance 
of the heat conducted to the droplet with the latent heat of the 
vaporization. The integration gives the following vaporization laws 
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These laws, which agree with those obtained by Tishkoff^, are valid for 
any value of £ , provided the bubble contains no droplet other than the test 

droplet. This standard law will be termed "canonical vaporization law," the 
spherical bubble surface will be identified by "canonical environment" and the 


bubble temperature by "canonical 

temperature. " 

In 

contrast to 

what 

is 

described above, an alternative 

law in which 

the 

vaporization 

rate 

is 


determined by the gas temperature on the surface of the transition sphere (see 
Fig. 2) will be called "renormalized vaporization law," the surface of 
transition sphere is "renormalized environment" and temperature on the surface 
is the "renormalized temperature" which is numerically equal to the local gas 
temperature of the continuum flow field. 
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Renormalized Droplet Law 


Since the continuum theory of sprays predicts the local gas temperature 

but not the canonical temperature, the canonical vaporization laws (5) and (6) 

cannot be used in spray calculation as previously described. The alternative 

proposed in this paper is to use a renormalized vaporization law. The 

renormalized law is described in this subsection. By adopting a mathematical 

procedure involving appropriate linear combinations of Eq. (2) governing 

and 9 one obtains two homogeneous equations governing 

m F T 

a + y and a - e • These two equations are integrated and joined with the 
T F F F 

inner solutions on the surface of the canonical bubble. The resulting 
solutions are given by 
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In order to obtain the renormalization laws, one first determines 

A A 

a and by replacing the v\ appearing in the lower limit of the 

Tco Fco 

integration of Eq. (7) by n * Subsequently, by substituting the resulting 

co 

expressions of and ct^ in Eq. (5) and (6) one arrives at the following 

r Tco Fco 

renormalized vaporization laws: 
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where C is the vaporization correction factor calculated as follows 
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The vaporization laws (8) and (9) remain incomplete until the distribution of 
y is provided. The determination of y requires the knowledge of the 
canonical temperature c* Tc (n) of a quasi-droplet located at n, see Eq. (12). A 
theoretical procedure of the determination of « Tc ( n ) tn terms of the local gas 
phase temperature, ci^n), is described in the following. 


Mean Canonical Bubble Temperature 
of a Field Droplet 

In the analytical estimate of a mean canonical temperature one assumes 
that each field droplet in the transition sphere of RND has a canonical 
structure, i.e., the temperature and concentration profiles are those given by 
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canonical model. Consider a canonical field droplet located at the field 
point, n , measured from the center of the test droplet. Let y be the radial 
location of an arbitrary point within the bubble of the field droplet measured 
from the center of the field droplet. Then, according to the result of the 
DIB model, the temperature distribution in the bubble of the canonical field 
droplet is given by 


a T (n,y) = y {expt—i — (1 - -)-l]} 
°(n ) y 


(14) 


where 

a = 4 7T p D r (n) 

Additionally, the canonical temperature, a T{; (n) is obtained by replacing 
y in Eq. (14) by y c ; i.e., the radius of the bubble. The result is 

a Tc(n) H a T (n ’ y c } = Vexp U -~>J - 1> (15) 

c 

Since the local gas temperature o T (n) equals to the volume averaged mean gas 
temperature in the bubble of a field droplet, one writes 

i y c 

a T (n) = — fj a T (n,y)y dy 

c (16) 
where V c is the non-dimensional volume of the bubble given by 
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v = 4 (y c 

c 3 c 

In order to express a Tc (n) in terms of a T (n), one inroduces the 
correlation function X as 


x(n,y c ) = a Tc (n,y c )/a T (n) 


07) 


By 


comparing Eqs. (16) and (17), one identifies A(n,y c ) as follows 
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Thus, the non-dimensional vaporization rate, p , of a field droplet is given by 
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where £ c (n)=£ co 0(n )g(n ) (20) 

By substituting (19) into Eq. (2), for 1 ^=^, one obtains the following 

equation. 
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( 22 ) 
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In the ensuing analysis, \ will be assumed to be unity. This corresponds 
to droplet vaporization with a high transfer number. The general case, 

A*l, can be solved by an iterative analysis with a guessed value of a(n) to 
calculate approximate <j) from Eq. (21). The iteration will continue until the 
vaporization rate prediced from the iterative solution converages to the 
guessed value of a. It is expected that the temperature rise in the radial 
direction is faster when \>1. Details of the general case will be reported in 
the future. 

The last step required for determining p which appears in the scaling law 

is to determine as a function of q. With A=l, one can show that the 

expression for <f> which satisfies 4> =4> co at n = n cQ and 4>=4' t . s at n=n ts is given 

by a linear combination of two homogeneous solutions W and W , that satisfy 

dw. dW~ 

the canonical boundary conditions; H,(n ) = 1, -z (q )=0 and -= (q ) = 1. 
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where 4> , the characteristic value, is determined in terms of <j> as 

co ‘•’ts 

follows. 

By equating the vaporization rate in the canonical form (5) and 
renormalized form (3) and by using the definition of ij> given by Eq. (22), one 
obtains 
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On substituting the expression of <t«(n) given by Eq. (23) into RHS of Eq. 
(25) and by inserting the resulting expression in the y-term appears in the 
denominator of the Eq. (24), one obtains 
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Finally, the scaling factor Eq. (10) can be expressed in terras of 


Wj and ^2 as follows 
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4 Numerical Results 

A numerical analysis will examine the structure and scaling laws of RND 

and their dependence on the principal collective interaction parameters of a 

stationary cloud of n-octane droplets with the following fuel properties; 

3 

p £ =707 Kg/m = 398.7 k, L = 71.7 Kcal/Kg, and W p = 114.14 Kg/kg-mole. 

(1) Pair-Distribution Function and Canonical Bubble 

In the absence of experimental data, a pair-distribution function is 
constructed on the basis of the geometrical distribution of molecules in a 
dense liquid. The following two parameter function is adopted for the 
numerical calculation. 


g(n) = 1 + a exp (-bn) cos (27rn co ), 2n co < n < n (29) 

where a and b are constants to be determined from the experimental data. In 
the present analysis a = 1.8, b = 0. 65 are chosen. The resulting distribution 
patterns are shown in Fig. 3. The corresponding signatures of the inverse of 

the radius of canonical bubbles with n cQ = 5, 10, and 15 are shown in Fig. 

4. 

(2) Vaporization Shape Factor 

A pronounced increase of y in the radial direction is observed for a 
smaller value of droplet spacing, i.e. = 7.5 in Fig. 5. The ratio of the 

vaporization rate of the droplet located at the first coordination shell, for 
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the case of n - 7.5, is approximately 5% of the corresponding value of RND 

with n = 15. This trend of a higher increasing rate of y at smaller 
co 

canonical bubble radii, or droplet spacing, is a common feature for small 
droplet spacing. This is confirmed for the cases characterized by pair- 
distribution functions with different values of a and b. The results suggest 
that "y -stratification" is a unique feature of RND in non-dilute clouds or 
spray. 

(3) Temperature Distribution 

High y- stratification at smaller droplet spacing, as illustrated in Fig. 
5 is due to the steep radial temperature gradients in the transition sphere, 
shown in Fig. (6). Indeed, the comparison of "y- T stratifications" suggest 
that (1) the rapid vaporization in outer layer of the cloud collectively 
quench the environment and thereby reduce the vaporization of the test droplet 
and (2) the increase in n ts at a fixed droplet spacing tends to reduce the 
inward heat transfer rate and thus suppresses the vaporization of the test 
droplet. 

(4) Vaporization Rate - Correction Factor 

The correction factors of RNDs for three selected values of n cQ , Fig. 7, 
are found to decrease monotonically as the group combustion number of RND 
increases. Saturation is projected to occur when G ^=30~40 with h cQ “ 7.5. 
While the group combustion number is a primary factor controlling the 
magnitude of C , Fig. 7 shows a small variation in the correction factor for 
two RND f s which have the same group combustion number but different 
renormalization number, 6 - n. /h • 

L S CO 
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Discussion 


Scaling Law with Linear Stratification fadel 

In order to facilitate the practical application of the scaling law, the cor- 
rection factor is integrated by adopting the following functional form of g u 


gP = gpK^n 


where is a stratification coefficient, and gp is the mean value of the 
weighted vaporization shape factor. The correction factor predicted by the 
linear stratification model described above, is given by 
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SPRAY CLASSIFICATION 

Numerical assessment of the scaling law suggests the following structural 
classification of non-dilute sprays* 

1. Diffusively Dense Cloud 

In a moderately dense cloud, RND is expected to have a large transition 
sphere that has no \i — stratification* The reduction in the vaporization is 
attributed to uniform thermal quenching* The correction factor given below is 
described by the group combustion alone: 

C v = (1 + V 1 (31) 
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2. Densely Stratified Cloud 


This cloud is featured with a strong stratification in a transition 

sphere that causes an intense collective quenching and a reduction in the 

vaporization rate of the test droplet* The renormalization number n /r\ is 

ts co 

larger than unity so that the correction factor is given by 


C 

v 



g r » (1 



ts 



(32) 


3. Sharply Dense Cloud - Fine Structure 

When the coordination shells contain the largest possible number of the 
droplet so that the renormalization number is not excessively large compared 
with unity, the correction factor depends on all the collective parameters: 
G™, n /n and n • Two sharply dense clouds with the same group combustion 
number will exhibit structural variation when the renormalization number is 
different. 

5 Conclusion 

The present theory portrays a droplet interacting with its neighbors by a 
minimum sized "dressed” droplet structured with a centrally located droplet 
and its bubble, e.g. Droplet-in-Bubble , surrounded by a quasi -droplet cloud 
spreading through a transition sphere. The distribution of quasi -droplets and 
the size of the cloud are described by a pair distribution function. The 

introduction of the transition sphere and the correspondence hypothesis 
postulating the equivalence between the average gas properties on the 
transition sphere with the gas properties of the continuum flow field, as the 
limit R tg / ^ goes to zero, constitutes the fundamental link between the 
discrete droplet dominated region with that of the surrounding continuum 
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flow. This feature provides the basis of the theory of short-range 

interaction and droplet rate processes. 

The vaporization scaling law reveals that the droplet vaporization is 
retarded primarily by the collective thermal quenching and the formation of 
high vapor concent rat ion in the transition sphere. The retardation is scaled 
by combined topological -thermochemical parameters: the group combustion 

number, renormalization number and droplet spacing. The scaling laws suggest 
that non-dilute clouds or sprays can be classified into: (1) "Diffusively 

Dense" clouds in which the vaporization reduction is attributed to the effects 
of uniform thermal quenching that can be scaled by the group combustion number 
alone; (2) "Densely Stratified" clouds, with marked radial stratification due 
to the large temperature gradient which exists in the transition sphere; and 
(3) "Sharply Dense" clouds, with a relatively large canonical bubble, and 
radial stratification so that the scaling law depends on G^, n t and B# 
Selected case studies show that the correction factor falls monotonically as 
the group combustion number increases and the saturation is expected to occur 
when the is of the order of 30 ~ 40, with the droplet spacing of 

approximately 7. 5 times of the droplet size and environment temperature of 
500 K. 

Finally, the validity of quasi-steady (Q-S) approximation is examined by 

comparing the relative magnitude of the characteristic time, t,.-_ , of 

dif f . 

difusion through the transition sphere and the life time, t-,- , of RND. By 

life 

adopting the renormalized vaporization laws, Eqs. (8) and (10), and a standard 
diffusion time for the transition sphere of RND, one concludes that the Q-S 
approximation is valid when the following inequality is satisfied 


90 


(33) 


t dif . 
t life 


2C n~ 1 (l + 
V ts n 


B) 


where B = C (T-T )/L. 

P b 


With p ~lo\g/ra'\ p~l kg/n^y T-T^~100°K L-lO^Kcal /kg , 

C -0.5 Kcal/Kg-K the numerical values of is 5 x 10 

p t 'lif£ 

Alternatively, by using the asymptotic form of for high G^, deduced 

from Eq. (30), i.e. , c ~G * = 4 G n gy ) , 1 one transforms Eq. (33) into the 

V KIM O S t S 

following form 


gyG > > G* (34) 

s s 

where G =12 In (1 + B) 

5 P £ 

3 

Recalling that G =bi\r\r n is equal to a third of the void fraction, and 
° s £ o 

gp~ 0(1), one concludes that the Q-S approximation is valid in high G-sprays 

* 

when the void fraction exceeds a third of the values of G • For typical values 

of p , p,T-T , C and L given above, Eq. (34) reduces to 

£ bp 

nr 3 > > 4xl0 -4 (33) 

io 

For example, with SMD=100y , the Q— S approximation is valid when n is greater 
than 3200 drops/cc, whereas for SMD=200y , the corresponding number density is 
400 drops/cc. 

The criterion (33) or its alternative form can be encoded in a spray 
numerical code to test the validity of Q-S assumption at each point of a spray 


91 


flow field. When the approximation is invalid, an alternative transient 
droplet laws should be used to determine the interfacial processes. Note that 
such droplet transient processes should be able to predict the change in the 
nature of interfacial processes. For example, a vaporizing droplet may reach 
the state of saturation when the droplet spacing falls below the critical 
value or when the thermal shielding retards the heat flux to the droplet. 
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ELEMENTAL DROPLET MODELS 
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Fig. 1 Droplet Models; Elemental Droplets 




COMPOSITE DROPLET MODELS 
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Fig. 1 Droplet Model; Composite Droplet 
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Funtion 



V 


Fig. 3 Pair-distribution function of quasi-droplets 
in transition sphere for three selected test 
droplet bubble sizes. 
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Fig. A Signature of the inverse of cannonical bubble radius 
for three selected test droplet bubble sizes. 
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Fig. 5 y-stratification; Vaporization shape factor distribution 
in transition sphere for three selected test droplet 
bubble sizes- 




Fig. 6 T- stratification = gas temperature distribution in 
transition sphere for three selected test droplet 
bubble sizes. 
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Fig. 7 Correction factor for vaporization rate vs group com- 
bustion number for three selected test droplet sizcMi 
Exact Solution. 
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Effect of Particle Velocity Fluctuations 
on the Inertia Coupling in Two-Phase Flow 

by 

Donald A. Drew 

Department of Mathematical Sciences 
Rensselaer Polytechnic Institute 
Troy, NY 12180-3590 

Introduction 

The fluid dynamics of flows of dispersed materials in a fluid is fundamental to suspen- 
sions, bubbly liquids, droplet flows, pneumatic transport, fluidization and erosion. Equa- 
tions of motion to describe these materials must deal with the interactions between them as 
well as the deformation of the carrier fluid. Models that treat assemblages of solid particles 
have been proposed and studied (Jenkins and Savage, 1983) that result m the particles 
behaving like a gas, with a pressure due to the fluctuations in the velocities that is at- 
tributed to collisional motions of the indivual particles. Models for fluid-particle mixtures 
(Drew 1983) do not include this effect (Passman 1989) The purpose of the present paper 
is to derive constitutive equations to supplement the equations of motion that include the 
effects of the particle velocity fluctuations. The particle motions are assumed to be at a 
sufficiently high Reynolds number that the fluid motion is inviscid, but viscous effects such 

as boundary layer separation are neglected. 

Equations of Motion 

The appropriate general average is the ensemble average. The ensemble average is 
the appropriate generalization of adding the values of the variable for each realization, and 
dividing by the number of observations. We shall refer to a “process” as the set of possible 
flows that could occur, given that the initial and boundary conditions are those appropriate 
to the physical situation that we wish to describe. We refer to a “realization’ of the flow 
as a possible motion that could have happened. Generally, we expect an infinite number 
of realizations of the flow, consisting of variations of position, attitude, and velocities of 

the discrete units and the fluid between them. 

If / is some field (i.e., a function of position x and time t) for some particular real- 
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ization, /<, of the process, then the average of / is 

(1) /(x,<)= f /(x,<;^)c!m(/i) 

Jm 

where dm(n) is the measure (probability) of observing process // and M is the set of all 
processes. We refer to M as the ensemble. The ensemble average is the fundamental 
average that allows the interpretation of the phenomena in terms of the repeatability of 
experiments. Any one exact experiment or realization will not be repeatable; however, any 
repetition of the experiment will lead to another member of the ensemble. 

In order to apply the averaging procedure to the equations of motion, we shall need 
some results about the averaging procedure. We shall also give these results for time- and 
volume averaging. 

In order to average to the exact equations, we need expressions for df /dt and V/. If 
f is ‘well behaved 1 ’, then it is clear from the definition of the ensemble average that 

,,, W = 3J 

dt dt 

and 


(3) 


V/ = V/ 


Functions are generally discontinuous at the interface in most multiphase flow. They 
are well behaved within each phase, however. Thus, consider XjtV/, where AT is the phase 
indicator function for phase k : 

(4) AT = ( *’ dxek ' 

1 0, otherwise. 

In the averaging process we will require the result 



dXk 

dt 


+ ■ VAT jt = 0. 


This relation has a reasonable physical explanation. Note that it is the “material” deriva- 
tive of A, following the interface. If we look at a point that is not on the interface, then 
either A* = 1 or A* = 0. In either case, the partial derivatives are both zero, and hence 
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the expression (5) is zero. If we consider a point on the interface, if we move with the 
interface, we see the function X k as a constant jump. Thus, its material derivative is zero. 

The averaged equations are 
Mass 

q Y 

( 6 ) + V • = p(v - Vi ) • VAT 


Momentum 

dXkpv 

( 0 XI ^ v • Xkpvv = V • X k T + x k pg + (pv(v - v t ) - T) ■ VAT-. 


Next, we define the appropriate average variables describing multiphase mechanics. 
First, the geometry of the exact, or microscopic situation is defined in terms of the 
phase function X k . The average of X k is the average fraction of the occurrences of phase' 
k at point x at time t. 

( § ) a k = X~k 


The quantity dX k /dn is the delta function defining the interface, its average is the inter- 
facial area density. 


( 9 ) 


a. 


dX k 

dn 


All the remaining variables are defined in terms of weighted averages. The main. 


or 


phasic 1 variables are either phasic weighted variables (weighted with the phase function 
A A;) or mass- weighted (or Favre) averaged (weighted by X */>). 

The “conserved” variables are the density 


( 1Q ) p x k =X k p/a k , 

and the velocity 


( 11 ) 


= X k py/a k pl 
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The averaged stress is defined by 


( 12 ) 


T* = X k T/a k 


The average body force is 

(13) ir = V^I. /a„pl 

As discussed above, several terms appear representing the actions of the convective 
and molecular fluxes at the interface. The convective flux terms are the mass generation 
rate 


(14) r * = p(v - v,-) • VAT 

and the interfacial momentum flux 

(15) vgr fc = pv(v - Vi ) • VAT 
The interfacial momentum source is defined by 

(16) M k = -T • VI, 

The motion of the interfaces gives rise to velocities that are not “laminar” in general. 
The velocity fluctuations may be due to turbulence or to the motion in the phases due 
to the motion of the interfaces. The effect of these velocity fluctuations, whatever their 
source, on a variable is accounted for by introducing its fluctuating field (denoted by the 
prime superscript), which is the difference between the complete field and the appropriate 
mean field. For example, 

v k = v - V 

Then, 


( 17 ) 


X k pvv =X k p(\ x k p + v'Xv*/ + v^) 
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The Reynolds stress is defined by 


(IS) Tf £ = X k pv' k v' k /a k , 

The averaged interfacial pressure , and shear stress t k% are introduced to separate 
mean field effects from local effects in the interfacial force. The interfacial pressure is 
defined by 

( 19 ) Pki = pdXk/dnk/a , 

and the interfacial shear stress is 

(20) T ki =TkdXk/dn k /ai 

Thus, 

M* = - T • VX k 

- r • 

=Pki VXT - Tki ■ VX~k - T', • VXk 

(21) =PtiVajt-Tjfc ) V«it + M' t) 

where we define the interfacial extra momentum source 

(22) M' fc = M k + p*,-V a k - Tjti ■ Va fc 


and introdxice 

+ r' fcl = -(p - pti ) I + (r - r*,)- 

The averaged equations governing each phase are 
Mass 


da t r t 


+ V • at -pl^Y 


Momentum 


dopin' 

dt 


+ V • a k ptn^Y 


= v • a k 


+ Tf e 


zPk g + M k + vj 
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The equation of conservation of mass for phase k (23) can be used in the momentum 
equation (24) to yield the Lagrangian form of the momentum equation: 

-s 

i r\ . — — — i n . 


a kPk 


Dt 


= OCkPk 


+v*/-Vv'/ 


) 


(25) 


= V-a^Tj+Tj^J + a^g + Mi 
+ P*iV a k - T ki ■ Va k + k: - v x k p )r k . 


The jump conditions are derived by multiplying the exact jump condition by nj • VA j 
and civeraging. This process yields the following conditions: 

Mass 


(26) 


hi + r 2 — o 


Momentum 


(27) M\ + M 2 + y^T 1 + v? t r 2 = m 

Constitutive Equations 

The exact equations of motion can be solved for the flow of an inviscid, incompressible 
irrotational fluid around an isolated sphere. We shall use this solution to derive information 
about constitutive equations for the force on the dispersed phase, the average stress, the 
Reynolds stress, and the interfacial pressure when the particle phase is allowed to have a 
random velocity. 

The fluid velocity at x for the irrotational flow of an incompressible inviscid fluid is 
expressed in terms of the velocity potential by 

(28) v(x) = V<£(x) 

The continuity equation becomes 

(28) 0- V v = V- V<f> = V 2 </>. 

The pressure at any point x is given in terms of the velocity by Bernoulli’s equation. 

^ / X 

(30) p- p(— + — | V^>| 2 ) = po = constant. 


108 



Consider a sphere located at a point z in a flow field, moving with velocity v p . The 
boundary condition at the surface of the sphere is 

(31) n • v p = n v = n at |x - z| = a, 

where a is the radius of the sphere, n is the normal to the surface of the sphere and v ?) is 
the velocity of the sphere. The boundary condition far from the sphere is 

(32) (f> -> as |x - z| — > oo, 
where 

. , , 1 

4 > oo = Vo(<) • X + -x • e f ■ X 

is the velocity potential that would exist in the fluid if the sphere were not present. Here 
v 0 (t) is the (unsteady) velocity of the fluid at the origin, and e/ is the rate of strain tensor 
for the fluid. We shall assume that ej is constant. 

A convenient form for the solution of this problem is given by Voinov (1973), and is 

, , , 1 

<P = V 0 (<) • x -| - -x • e/ ■ x 

+ \ ( v p - v o(0 - z • e f ) ■ (x - z) j 

(33) + I( x _ z ). e/ .( x -z)(^ 

If there are many spheres in the flow field, the solution given above (33) will still be 
a good approximation if the spheres are sufficiently far apart that the flow disturbances 
due to the individual spheres do not interact. That is, the flow must be sufficiently dilute. 
Thus, we can think of each sphere as “isolated” in the sense that it only interacts with its 
neighbors through the averaged fields. We assume that each sphere lies in a “cell,” and 
inside that cell, the velocity is given by eq. (33). We shall approximate the cell to be a 
sphere of radius R. We choose R so that 

^7ra 3 /^rf = a. 
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Averaging 

We now introduce a means of evaluating the averaging process. The first aspect of 
the ensemble average is that the sphere velocity is random, with a distribution function 
f 1 1 )( Vp, x, t ). The second aspect is that the sphere and the surrounding cell can lie with 
the sphere center anywhere within radius R of the space point x. The average is performed 
by first integrating over the distribution of the velocities that the sphere can have, followed 
by an integration over the possible positions that the sphere center can have. Note that if 
|x — z| < a, the material making up the sphere occupies the field point x and if |x — z| > a 
the fluid occupies the field point x. 

The average of a quantity g(xj;z,v p ) is performed in two parts, that is first, we 
perform a conditional average of g for given sphere position z, integrating over the velocity 
space v p , followed by the average over the spatial positions the sphere can have. We 
assume that the distribution of positions is such that 

dV r , Vofd(x,<) ‘ 

§7T if* _ X Otd(x,t) 

is the probability of finding the sphere in a volume dV surrounding the point z, where 
x* — x — z. Thus, for the average over the fluid of a quantity < 7 , we have 


(34) 


g(x,t\z) = / g(x,t; z, v p )/ ( 1 ) (v p , z, t)dV v 

Jr* 


Here the notation g{x,t |z) is intended to suggest the conditional average assuming the 
sphere is located at z. The average of g over the fluid phase is then given by 


(35) 




1 7 r(f ? 3 — a 3 ) 


f f f g(x,t\z)dSldr, 
Ja J Jn(r) 


where Q( r ) is the sphere of radius r centered at x, and the integration is over the z variable. 

It will be convenient to introduce the average particle velocity and the fluctuation 
particle velocity as 


(3G) 


Vp(x,#)= f v ;) / ( 1 ) (vp,x,f) dV Vp 
Jr 3 


(37) 


Vp(x,t) = Vp - Vp(x,<) 


no 



Note that 


(38) \' p = 0 

The particle kinetic energy per unit particle mass is 

( 39 ) u d e (*, t ) = \j |vj,| 2 / ( 1 ) (v p , x ,f) dV v 

1 Jr 3 

and the Reynolds stress for the particles is defined by 

(«) Tj*(x,<)= -Pd f vX/l 1 >(v p ,x,()rfV' v , 

JR 3 

In order to evaluate the integrals appearing in the averaging process, we must express 
the z dependence of the velocities in terms of x and x' = x — z. We have 

V/(z) = vy(x) -x' • ey 

and 

Vp(z) = Vp(x) - x -e p 

where is the velocity gradient tensor for the average particle motion. We shall assume 
that this tensor is constant and symmetric. 

We have 


v(x, t; z, Vp) = V<£(x, t: z, v p ) 

= v /( x ) + J( v /( x ) “ v p (x) - vj, - x' • (ey - e p )) ^ 

- ^( v /( x ) - Vp(x) - Vp - x' • (ey - e p )) • x' x' 


(41) 


2 , /a 5 

+ 3 xe 'b 


S , 

-x • er ■ x 
3 ' 


a 

k7 


Note that vy(x) is the fluid velocity that would exist at x if the sphere were not present, 
and v y(z) — v p is the relative velocity between the sphere and the fluid evaluated at the 
sphere center. It is convenient to have expressions for the integrals of powers of x' over 
Of /• ). For these integrals, we note that 


(42a) 



. . x'd£l = 0 
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if the factor x' appears on odd number of times, and 


(42b) 

(42c) 

(42d) 

where E is a fourth order isotropic tensor defined in Cartesian coordinates by 

'Eijkl = fiijfikl + + fiilfijk- 

We further note that if v is a vector, and e is a symmetric second order tensor with e,; = 0, 
then 

Ej j kiv jf^ui — 2vjCji- 

Derivation of Averaged Quantities 

In order to average eq. (41), we note that the average over the velocity fluctuations 
gives no contribution, by eq. (38). Then averaging over z gives 

(43) v^(x,i) = v/(x,t). 

Similarly, to obtain the interfacial averaged velocity of the fluid, again the integration 
over the velocity fluctuations gives no contribution, and we have 



Substituting and performing the integrations lead to the result 
(44) v ci (x,i) = v/(x,f). 

This result is a little surprising at first. The fluid at the surface of the sphere satisfies the 
condition n • v = n • v p , but is allowed to slip in the tangential direction. After the passage 


/ dQ = 
Jn(r) 


= 47r r 


In(r) 


x'x'di} = -7rr 4 I 

3 


Jn(r] 


x'x'x'x'dQ, = — 7rr 6 E 

15 
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of the sphere, the fluid that was momentarily in contact with the surface of the sphere is 
again moving with the fluid. The result says that even during the time that it is in contact 
with the surface of the sphere, its average velocity is still equal to the average velocity of 
the fluid, and not of the sphere. 

Now let us compute averaged pressures using this formalism. The exact pressure can 
be computed by Bernoulli’s equation (30). In order to evaluate the derivatives in eq. (30), 
we note that x is constant during t derivatives, but dz/dt = v p = v p (z, t) + v p . Also, 
when evaluating V</>, both t and z are held constant. The pressure is given by 

(dv 0 1 [dvo 5v P (x,t) 

(45) p(x,i;z,v p ) = p 0 -Pc l ~Qf x+ 2 dt dt 


+ (v p (x, t) + Vp) • e f - (v p (x, t ) + v'p) • e p 



■^( v /(x) — (v p (z,f) + Vp) — x' • (ej - e p )) • (v p (z,f) + v p ) ( 3 


-|(v/(x) - (V p (z, t ) + vp - x' ■ (e/ - e p )) • x' 



x' • (v p (z,<) + v'p) 


-J( v p(M) + v p) • e / • x ' ’ e f ' x ’ x ' ( v p( z ’ < ) + v p) (r 5 ) + 2 V/ V/ 



+ ^l v /( x ) ~ ( v p( z ^) + v p) — x/ ' ( e / e p)l 

o 


+^ v / • ( v /( x ) - ( v p( x ’ 0 + v p) - x ‘ ( e / - e p )) ( ^3 

+ — ((v/(x) — (Vp(x, t ) + Vp) — x • (ej — e p )) • x ) 

8 



10 


+^( v p( x ^) + v p) ■ e / ■ x 


1 


- |(v p (x,f) + v' p )-x'x' - e/ x' 




g(v/(x) - (Vp(xrf) + vj,)) • e/ • x' ( ^ 


-|(v/(x) - (v p (x, t) + vp) • x'x' • e f - x' 


+ |( v /(x) - (v p (x,t) + vp) - x'x' -e/ -x' ^ 
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where we have ignored terms of order c ^ , e p , and eye p . Averaging over the velocity 
fluctuations gives 


-f .1 X l'dVQ 1 

P(x, t|z) = Po -Pc -=7- -X + - 

at 2 


dv 0 dv p (x, f) 
dt dt 




+v^(x, t.) ■ ej — v p (x, #) • e p ] • x' 


q(v/(x) - v p (z, <) - x' • ( e/ - e p )) • Vp(z, t) - uj? (z, f) 


|(v/(x) - Vp(z,0 - x* • ( e/ - e p )) • xV-v p (z,t) + ' x ' 

2 p d 

-\v P {z, t) - e/ ■ x ; + ^x' • e f • x'x' • v p (z, t) ( 4 ) + ~v f ■ 



-v/ • V/ 


+ 


§M X ) ~ v r( z , 0 - x' • (e, - e p )| + ^f(z,t) 
+ J V / • (v/(x) - v p (x,i) - x' • (e/ - e p )) ( ~ 





+ 


((V/(xj 


v„ X 


,0-x'.( e/ -e p ))-xT- 


,\2 9 x'-T«‘(z,t). x ' 


Pd 


a 6 

TTo 


+ 3V p (x,t)-e / -x' j - jv p (x,t,) - x'x' -e f -x' [ ^ 

^(v/(x) - v p (x,t)) • e 7 x' 

-^(v/(x) - v p (x, *)) • x'x' • e f • x' 

+ ^( V /( X ) - Vp(x, <)) • x'x' • e/ • x' ^ "-j j 

The spatial integration is tedious, 1> 1 1 1 results in 




(46) 


dx i 


Pc - Po - Pc-Qf • X - 9 V /(x) • V/(x) 
where we ignore terms of order a/ R in addition to those ignored previously. We also obtai 
( 47 > Pci(x,i) = Pc - jp c |v c (x,f) - v d (x,0| 2 + \p c u% e . 


■mi 
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The interfacial force density M d is given by M rf = pV X d . Thus, 

1 r / \ 

(48) M 4(x,0 = (1 - — J dSI 

This can be computed by substituting eq. (45) for the pressure, and recognizing that 
n = x'/a. We must also expand the terms u^ e (z,t) = u^ e (x,i) - x' • Vu^ £ (x,t) and 
T‘j e {z,t) = T £ e (x,f) - x' ■ VTf(x,/) The result is that 


(49) 


+ &dPc ( ey 


Mj(x,t) = p z c ^ad 

dv^(x,t) _ dvYix rt + ^xp (x> f } . Vvgp(x, /) - v7(x, t) • Vv5 P (x t <) 


dt 


dt 


-l(v^(x,<) - v7(x,/)) • (Vv^(x,<) - Vv*/(x,0)) 


+p* ( i(v? p (x, t) - v7(x, *)) • (v*'(x, 0 - v-(x, t))Va d 


T^v^x,*) - v 7 (x,<)(v^(x ,0 - v'/(x,0) • Va ^) 


-atflzVv? - • Tf - ’ T 


-X ' „.fie 


P? 9 


I'dP, 


20 


Pi 20 


20 


Pd 20 


^/?e 

cJ 


Note that no drag force is present in eq. (49). This is the result of D Alembeit s paradox, 
that is, there is no net force on a body moving at a constant velocity through an mviscid 
fluid at rest. 

If a distribution of stresses is applied to the surface of an elastic body, there results a 
distribution of stresses inside the body. These induced stresses are important in computing 
constitutive equations for solid-fluid mixtures. The average stress inside the sphere is given 


by 


T d {x,t) = 


7 rcr 


rn 

J0 J Jd(r) 


T(x,f|z)dO 


where 

T(x,t|z) = / T(x,f;z, Vp)/ ( 1 ) (v ?) ,z,/) d\ Vp 

Jr 3 

Here T(x,t; z, v p ) is the stress at point x inside a sphere having its center at z at time t. 
We shall assume that the spheres are linearly elastic solids, but we shall assume that the 
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deformation is sufficiently small that the fluid motion is unaffected by the deformation of 
the spheres. Then the stress-strain relation is given by 

(50) T = p s [Vu + (Vu) lr J -f A S V • ul 

This can be written as 


(51) 


T = <7 + 01 


w 


here 


© = ( + -p s ) V - u 


° - Us |[Vu + (Vu) ir ) - • ul 

The spherical part of the stress satisfies (Love, 1932) 


(52a) 


V 2 0(x, t] z, v p ) = 0 


( 52b ) 0(x, t; z, v p ) = -p(x, t; z, v p ) on |x - z| = a 

Averaging over v p gives 

(53a) V 2 0(x, f |z) = 0 


(53b) 0(x, t\z) = — p(x, t\z) on |x — z| = a 

Solving and performing the integration over z gives 

(54) @d(x,t) = -p ci (x,t). 


The solution for a is also given in Love (1932) and can be averaged in v p and then 
integrated in z to give 


*5(x,<) = P 


~ go ( (v-(x, t) - v7(x, t )){y* e ' (x, t) - v7(x, 0) 


I ?\ 

Pd ) 


(55) 


+ ^ (K P (x, t) - v x /(x, t ) | 2 + 2 u* e ) I 
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We next turn to computations of the Reynolds stress using the velocity fluctuations 
due to the inviscid flow around a sphere. Using the expression for the velocity (41) and 
the average fluid velocity (43), we see that 


(56) 


v' c (x,t;z,v p ) = ^(v f ( z) - v p ) 



- |(v/M - Vp) • x' 



+ 3 X 6/ 


' / t 

x -e r x 



x 


/ 


Averaging over the particle velocity fluctuations yields 


Tf e (x,f|z) = -p c Vc(x,t|z)v' c (x,f|z) 


+ Pc 


1 {a*\ Tf 


Pd 


4 \ r 




t He 




-) + (*' 


rp Re 

— )*'] 
Pd 


9 / a 6 \ , , , Tf 1 

+ - — x'x' x' • * x ) 

4 \ r 1 0 / Pi J 


The integration over z can be performed, yielding 

1 


Tf(x,<) - 


— X 
C 


nn He 

(v^-v7)(v^-v7)-^- 

Pd 


(57; 


+3((vr-v7)-(vc p -v7) + 2uf)l]. 


The fluid fluctuation kinetic energy is uf e = • v' r , and can be computed by taking 

the trace of eq. (57) for Tf e . The result is 


(58) 


- v 7 1 2 + \ a * u d l e - 


Conservation of Fluctuation Kinetic Energy 

In analogy with statistical mechanics for assemblages of molecules, the theory of av- 
eraging as applied to multiphase flows allows the computation of averaged equations foi 
higher moments of velocity and pressure correlations. 

We start with the derivation of averaged equations for the fluctuation kinetic enei gv 
for each phase. The exact momentum equation is 

(59) + V • pvv = V • T + pg 
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We shall derive an equation for the evolution of the kinetic energy. The Lagrangian 
form of the momentum equation is 


(GO) 


= V • T + pg 


If we take the dot product of v with eq. (60), we have 


(61) 

We note that 


i „,2 


djv 

dt 


+ v • X-v 2 = v • (V ■ T) + pv • g 


v • (V • T) = V ■ (T • v) — T : Vv 
If we also return to the Eulerian form, we have 

1 


(62) 


d\pv 2 j. ,, 

— ^ h V • -pv 2 v = V ■ (T • v) — T : Vv + pv ■ g 


If we apply the ensemble average to eq. (62), we have 


dXkhpv 2 „ 

2 TV • X k -pv 2 \ = V • X k T • v - ATT : Vv 


dt 


1 

2 ' 


(63) 


+ X k pv ■ g - [p-v 2 (\ - v.) + T • v) • vat; 


We define the fluctuation velocity of phase k by 


(64) 


Then 


(65) 


v' =V-V- 


u 2 = {v' k f + 2v' k • vf + {vYf 


so that, noting that X k p\' k = X k p\ — X k p\ x k p = 0, we have 


Xk p\v 2 = AT p(v' k ) 2 + X k p2\' k • v x k p + AT p(v xp ) 2 = a k p k u^ e + a h p x k \ v xp 
Furthermore, 


t- 3 v = K)v t + + 2vt ■ vj'vp + avi ■ vj'vj + + (vvfn 


TT X <> 
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so that 


=x k p'- (v’ t fV k + + v ‘t ■ *« w l v * + a '^5«") 2v 




=a k q k + ®kPk u k ev k P - a * v k P ‘ T fc e + Q kPk 0 ( v k P ) 2 v k / 


Also, note that 


T v = T vl p + = T • v',. 


Note further that 


T : Vv = T : Vv ? \ p + T : Vv' fc 


Then we have 


X k T • v = X k T • v x k p + X k T • v' fc = a k T T k • v* p - a k ql 


and 


ATT : Vv = ATT : Vv* p + X k T : Vv^. = a k T x k : Vv 1 / + D k 


where = q p . + q£ and D k — A^T : Vv/ 
Next, in. the interfacial terms, we have 


>— u 2 (v - Vi) • VAT = htr p ) 2 n + v* p • v/rv + /?^K) 2 (v - Vi) • VAT 


and 


(T • v) • VA k = (T • v* p ) • VAT + (T • v/) • VA* = M k ■ v* p + W k 

The equation for the conservation of fluctuation kinetic energy then becomes 


d_ 

dt 


«kPl ( "F e + 


) +V • a k pl t + ^(^D 2 ) v^-V- (a k v T k p ■ T[ le ' 


- V • a k ( q/ + qj[.') + V • o,T^ ■ v/' + a^v/' • g 

- a k % : Vv 1 / - D k + |(^ p ) 2 r ft + v£ P ’ 

+ p\{v' k ) 2 {v - v «) • va t + M * • + Wk 




(CO) 


119 


The equation for the average kinetic energy is 


d^nj(vVi 

dt 


+ v ^ wr -(^) 2 = 


(07) 


v^- V-a fc (T fe +Tf e ) + M fc -v* p 

+ <*kpkg-v x k p + v x k p ■ v^k. 


Subtracting this from eq. (66), we have 


v ■ c k r k u^vY = ai Tf' : Vv'' - V • a t (qf + qj) 


i Re 


dt 


( 68 ) 


+ ^ V ) r * - P2 (u *) 2 ( v " v *') ' VX * + ^ ^ 


This equation has some interesting interpretations. First, note that the dissipation 
due to the Reynolds stress : Vv^T acts as a source of fluctuation kinetic energy, 

while its counterpart for the molecular dissipation : Vv^' does not appear in this 

equation. Dissipation on the macroscopic scale, then, winds up as different things on the 
microscopic scale. Also, the dissipation due to microscopic velocity fluctuations Dk implies 
a loss of fluctuation kinetic energy. Thus, loss mechanisms, such as inelastic collisions or 
viscous dissipation in the velocity fulctuations, cause a loss of fluctuation kinetic energy 
to heat. Finally, the working of the fluctuations at the interface, Wk appears as a source 
of fluctuation kinetic energy. 

Since this equation is unnecessary for the fluid phase, we shall ignore it for k = c. 
For k = d, we note that D d = 0 is consistent with the linear elasticity assumption and the 
assumption that the particle radius a does not change. Furthermore, if assume no phase 
change (T*. =0), and we ignore triple correlations in the particle velocity fluctuations 
( v ', v pVp = 0), we have 

du Re 

m • Vufvy = c d Tf : Vv7 + W d 


Discussion of the Force on a Sphere 

The equations of motion for the mixture are 


(70) 


da d 

dt 


+ V • a d v r d p = 0 
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(71) 


dot c 

~dt 


+ V • a c v xp = 0 


(72) a d ~p x d {^+^7 - W c") =-<*dV(p x c + \p>d*-\p X c |Vc P - v 7l 2 ) + V ‘ a <iT 


Re 

d 


+ v • OLd { P\ 


~ x 

c 


(v c (x,i) - Vd(x,t))(v c (x,t) - v d (x,0) 


+ — K p (x,0- Vj p (x,<)| 21 


+«#; Q + 0 ' vv '' ,(x ' () ~ v * f(x ' t} ' vv; " (x ’ 0 

-- [v; p (x, t ) - v'/(x, 0] • [Vv7(x, () - Vv7(x, 0l) 


+P x c (J(v* p (x, *) - v7(x, 0) • (v-(x, <) - v7(x, t))Va d 
+ ^(v7(x.<) - v7(x, 0)(v7(x, t) - v7(x,0) ■ Vaj) 


+/5?^V(a JU f) 


(73) 




1 


f dyg_ 

dt 


+ V ' { ~OQ adf> C 


+ V 1 / ■ Vv* p = -a c Vp; 


(v^-v7)(v^-v7) 


rpRe 

- 1 - d 

Pd 


-<*dpc I 9 


+3 ((v7 - v7) . (v7 - v7) + 2uf ) i] ) 

+ (i£|v7(x,0 + v7(x,f)| 2 + ?') 

( I [ <?5!£M) ^2 x,t) + v7(x, o • Vv7(x, i) - v7(x,() - Vv7(x, () 

\2 at 

[v7(x,i) - v7(x,()] . (Vv7(x,l) - Vv7(x,()l) 

-p’ fg(v7(x, () - v7(x, *)) ■ (V7(x, () - v‘/(x,t))Va d 
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+ ^(vj'(x,i) - v7(x, i))(v7(x,() - v7(x,0) ■ Vcfj 

+K^V(a«,«' ! ') + ■ («rfTj c ) 

It is also possible to calculate the force on a sphere at z by computing 



where the integration is over the variable x', with x = z + x'. This results in 


(74) 


Fp(z) = f-i ra 3 pf 


<9v y 

~di 


+ V/ . e/ + - 


9vr <9v,. 


at 


ot 


+ v, • e f 


Note that this force agrees with Taylor’s (1928) calculation of the force necessary to hold 
a sphere at rest in an accelerating stream, obtained by setting = 0 and v p = 0. The 
force is 


(75) 


§ (|™ 3 ) Kn- Vv / 


If we first take the gradients involved in eq. (72), using V • v)/’ = ( 1 / cv ,i ) ( <9a,/ / dt + v 1 ■ 
Vo rf ) and V -v* p = -( l/oc^do^/df + \ x c p ■ Vaj) , then set \ x f = 0, u£ c = 0, a d = const., 
and 'I'//'' = 0; and assume one-dimensional, steady flow, then eq. (72) reduces to eq. (75) 
in the limit as nj — > 0. Moreover, it is clear that it should. Consider the one-dimensional 
situation pictured in Fig. 1. The continuum model for the particles between x and x + Ax 
gives the rate of change of the momentum of the particles and parts of particles between 
x and x + Ax, denoted by Pd(x, x + Ax), as the stress force transmitted to the particles 
by the particle parts outside of the interval, denoted by (a^i - T^) | x + (ctrf( — i) • T ( ) ) | x -(.^ r , 
plus the force transmitted to the particles through their interface, denoted by MjAx. 
Thus, 


(76) Pd(x, x + A x ) = (a«/i • Tj) | x + («,/(— i) ■ T;)) | x +Ar + MjAx 


The sum of the forces on all particles with their centers in the interval from x to 
.r + A.r is equal to the sum of the pressure forces on all the particles involved. This is 
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denoted by S J pndS. This is equal to the rate of change of the momentum of all the 
particles, denoted by P p . Thus, 

pndS 

We note that eqs. (76) and (77) differ in the way they treat the particles being cut by the 
surfaces at x and x + Ax. The relation is that the 
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is the sum of the pressure forces on the surfaces of the cut particles at x whose centers arc 
outside the interval from x to x + A:r. A similar interpretation is valid for the cut particles 
at x + Ax. 

The terms on the right hand side of eq. (78) represent the resultants of forces on cut 
particles. If the approximate equation of motion inside the cut particle is V ■ T = 0, then 
the pressure force over the curved side, plus the stress resultant force over the flat side 
must add up to 0. (Note that if the particles are accelerating, then the forces add up to 
be the volume of the part of the cut particle, times the acceleration of its center of mass. 
Pr esumably, this force is small.) 

Conclusion 

Consistent forms for the interfacial force, the interfacial pressure, the Reynolds stresses 
and the particle stress have been derived for the inviscid, irrotational incompressible flow 
of fluid in a dilute suspension of spheres. The particles are assumed to have a velocity 
distribution, giving rise to an effective pressure and stress in the particle phase. The 
velocity fluctuations also contribute in the fluid Reynolds stress and in the (elastic) stress 
field inside the spheres. The relation of these constitutive equations to the force on an 
individual sphere is discussed. 
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ABSTRACT 

A model of dilute gas-solid flow in vertical risers is developed in which the particle 
phase is treated as a granular material, the balance equations for rapid granular flow are 
modified to incorporate the drag force from the gas, and boundary conditions, based on 
collisional exchanges of momentum and energy at the wall, are employed. In this model, it 
is assumed that the particle fluctuations are determined by inter-particle collisions only and 
that the turbulence of the gas is unaffected by the presence of the particles. The model is 
developed in the context of, but not limited to, steady, fully developed flow. A numerical 
solution of the resulting governing equations provides concentration profiles generally ob- 
served in dilute pneumatic flow, velocity profiles in good agreement with the measurements 
of Tsuji, et al . (1984), and an explanation for the enhancement of turbulence that they ob- 
served. 


INTRODUCTION 

Gas-solid flows satisfy the principles of mass and momentum conservation. While 
the Navier-Stokes equations govern the motion of the gas phase, there is not yet unani- 
mous agreement upon the form of the equations for the particle phase. Treating this phase 
as a continuum provides a framework in which techniques such as volume averaging may 
be used. This two-fluid approach results in the derivation of partial differential equations of 
motion. Depending on the situation, the stress tensor for the particle phase can then be 
modeled in order to close the system. 

For small particles, dilute in a turbulent gas, momentum transfer in the particle 
phase is due to turbulent diffusion of the particles. In this regime, Elghobashi & Abou- 
Arab (1983) have rigorously derived a two-fluid k-e model that relates the Reynolds stress 
of the particles to gradients of the mean particle velocity using an eddy viscosity that is a 
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fraction of the eddy viscosity of the gas. This model predicts the reduction of turbulent en- 
ergy observed in the presence of small particles (e.g., Modaress, et al. , 1984). Other mod- 
els for this regime include the works of Pourahmadi & Humphrey (1983), Chen & Wood 
(1985) and Berker & Tulig (1986). However, these models do not apply to gas-solid flows 
with large and heavy particles. 

For large particles, the experiments of Soo, et al. (1960) have shown that the in- 
tensity of particle velocity fluctuations may exceed that of the fluid, an observation that 
cannot be explained by treating the particle response to turbulence. Min (1967) attributed 
this high particle "turbulence" to particle-wall collisions. In this context, Lourenco, et al . 
(1983) have successfully modeled the air flow in a 10cm wide horizontal duct loaded with 
500)J.m glass particles. By analogy with molecular dynamics, these authors treat particle 
collisions using a particle velocity distribution function that satisfies the Boltzmann trans- 
port equation. They assume that the distribution is determined by particle collisions rather 
than by the gas turbulence. In other words, the gas affects the mean velocity of the particles 
but not their random motion. This assumption is justified in Lourenco's model, because the 
ratio of particle relaxation time to a typical large eddy turbulent timescale is large (Hinze, 
1972). Clearly, the success of Lourenco, et al . encourages further studies of the regime 
dominated by particle collisions. 

Another shortcoming of the two-fluid approach is the lack of attention given to the 
formulation of correct boundary conditions for the particle phase. Unlike the fluid phase, 
particles can slip at a boundary. For dilute suspensions of small particles, Soo (1969) sug- 
gested, without derivation, a set of boundary conditions by analogy with rarefied-gas dy- 
namics. Clearly, for regimes where collisions dominate momentum transfer among the 
particles, the boundary conditions of the particle phase should be carefully considered. 

In this study, a two-fluid model of the vertical flow in a gas-solid riser is proposed. 
For simplicity, the model assumes that the flow is fully-developed and steady. Particles are 
assumed to be sufficiently large or heavy so that momentum transfer in the particle phase is 
due to particle collisions alone. The granular theory used to describe the particles is related 
to the kinetic theory of gases, but it is not limited to dilute situations or to elastic particles 
(e.g., Jenkins, 1987). The conservation laws for mass and momentum in the particle phase 
have familiar forms. An additional balance law governs the measure w = V<v' 2 >/3 of the 
velocity fluctuations v’. The quantity w is the analog of the molecular temperature. For the 
particle phase, constitutive relations for the stress tensor and the flux of fluctuation energy 
are those supplied by the kinetic theory (Chapman & Cowling, 1970) and an existing form 
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for the collisions dissipation of fluctuation energy is employed (Jenkins & Savage, 1983). 
Finally, boundary conditions for the mean particle velocity and the granular temperature at a 
solid surface follow from a detailed description of the particle collision dynamics at the 
surface (Jenkins, 1988). Such rigorous derivation of the constitutive relations and bound- 
ary conditions are a major advantage of the granular flow theory against ad hoc continuum 
models for the particle phase. 

In this paper, the granular flow equations given, for example, in Jenkins (1987) are 
modified to include contributions from the gas. In particular, a drag term is added to the 
momentum equation and viscous dissipation of the particle rms fluctuating velocity is 
introduced. The resulting equations, constitutive relations and boundary conditions are 
presented for a dilute, steady, fully developed, upward flow in a vertical pipe with a circu- 
lar cross-section. A numerical solution of these equations is obtained and compared with 
the experimental data of Tsuji, £tak (1984). 


THF. TWO-FLUID MODEL 
Balance Laws 

The equations of motion for a steady, fully developed, axi-symmetric upward flow 
of a dilute gas-solid suspension are: 


Gas momentum 

Ta?< rt) - e a£' epg ' F=0 ’ (1) 

where e is the volume fraction of the gas, r and z are radial and vertical coordinates, X is the 
mean gas shear stress, p is the gas pressure, and F is the force per unit volume exerted by 
the gas on the particles; 


Particle momentum 

IA(rS)-(l-£)Psg-a-e)gf+F = 0 


( 2 ) 


and 

dN-0 O) 

dr _u ’ 

where S is the particle shear stress, p s is the density of the particle, g is the gravitational 
acceleration, and N is the particle pressure; 
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Particle fluctuation energy 


r"3r ( r q) + ^ cF ~ D - D' = 0 , 


(4) 


where q is the flux of fluctuation energy, v is the mean particle velocity, and D and D' are, 
respectively, the rates of collisional and viscous dissipation per unit volume. 

Constitutive Relations 

1) Volume supplies 

Equation (1) is the vertical component of the momentum equation for an incom- 
pressible gas exerting drag F on the particle phase. The mean drag is equal to the drag force 
on a single sphere based on the mean velocities, multiplied by the number of particles per 
unit volume: 


F = Cd I u-v I (u-v) (1-e) f(e) , 


(5) 


where u is the mean interstitial gas velocity, d is the particle diameter, and the drag coeffi- 
cient Cd is given in terms of the particle Reynolds number. Rep by 

Q = ^[. + 0..5RC»^ 7 ], ' (6) 

with 


Rep — 


lu-vlpd 


(7) 


in which p and p are, respectively, the density and viscosity of the gas. The drag coeffi- 
cient Cd is taken from Boothroyd (1971) and it is valid for 0 < Re p < 1000. The function 
f(e) is an empirical correction to the drag force on a single particle that allows for the pres- 
ence of other panicles. According to Foscolo and Gibilaro (1984), for the same velocity 

difference, the drag force per panicle increases in the presence of other particles according 
to 


f(e) « e 2 e- 3 - 8 , 

where the £ 2 term accounts for the ratio of superficial to interstitial velocities. 

Equations (2) and (3) are the vertical and radial components of the balance of parti- 
cle momentum in this simple flow. The shear stress and particle pressure result from trans- 
port of momentum between collisions exactly as in a dilute gas. The vertical balance in- 
cludes forces due to gravity, gas pressure gradient, and drag. 
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The first term of equation (4) is the diffusion of particle fluctuation energy. The 
second term is the energy produced by mean shear, the third is the rate of energy dissipa- 
tion per unit volume due to inelastic collisions. In the dilute limit of granular flow, the third 

term is given by 

D _ 240-e)£s w 3 (i_ e) 2 1 (9) 

yn d 

where e is the coefficient of restitution for a particle-particle collision; e=l for elastic parti- 
cles. This is obtained by considering the kinetic energy lost in each collision, then averag- 
ing over all possible collisions (Jenkins & Savage, 1983). 

The fourth term, D', is the rate of energy dissipation per unit volume arising from 
the drag force on the fluctuating particles. To calculate it, we first write the rate at which 
energy is gained by a particle as the scalar product of its total velocity and the particle drag 
force. Ignoring the velocity fluctuations in the gas, this is 

C d I u - v - C I Ijj-(l-e) ( U - v - C )•( v + C ) , 0°) 

where C is the velocity fluctuation of the particle. Then, upon evaluating the second term in 
the brackets in equation (6) at the mean velocity difference and averaging over all possible 
C, we find that the rate of energy loss per unit volume from the fluctuating motion is 

D’ = ^jCdl u-v I (1-e) w 2 . 


2) Fluxes 

In most situations of practical interest, the flow Reynolds number is high enough 
for the gas to be turbulent. For simplicity, the mean gas shear stress is modeled using an 

eddy viscosity Pt> 

. . du (12) 

x = e (p + (it) dF • K 

In this work, (i t is approximated using a polynomial fit to the experimental results for tur- 
bulent gas flow in a pipe (Hinze, 1975, Sec. 7-13). Thus, we assume that the gas 
Reynolds stress is unaffected by the presence of particles. We note, however, that several 
experimental observations have shown that large particles enhance the turbulence, even for 
conditions as dilute as (1-e) - 0.1% (e.g., Tsuji, etal, 1984 and Lee & Durst, 1982). Our 
view is that any treatment of the gas turbulence should not ignore the details of the parti- 
cles’ motion. In particular, the "pure" turbulent fluctuations in the gas must be distin- 
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guished from the randomness induced by particle fluctuations and the resulting change of 
the structure of turbulence must be quantified. Here, for simplicity, we ignore the influence 
of these induced fluctuations on the eddy viscosity. For a discussion of these effects, see 
the experiments of Boothroyd (1967), Nouri, gtal. (1987), Tsuji, etal (1984), and 
Modaress, et al. (1984), and the theoretical discussions of Owen (1969), Elghobashi & 
Abou-Arab (1983) and Genchev & Karpuzov (1980). 


The constitutive relations for the particle shear stress S, the particle pressure N, and 
the flux of fluctuation energy q for dilute flows of smooth, nearly elastic spheres, are iden- 
tical to those provided by Chapman & Cowling (1970, Secs. 7.4, 7.41, and 10.21): 


c S'frc , dv 
a = — Psdw ^r. 


96 


(13) 


N = (1-e) p s w 2 , 


(14) 


and 


q 


25VtT 

64 


p s d w 2 


dw 
dr ' 


(15) 


The shear stress is Newtonian and the particle viscosity is proportional to the rms fluctuat- 
ing velocity and the particle diameter. The energy flux is proportional to the gradient of 
3w 2 /2, the fluctuation energy per unit mass, and the corresponding coefficient depends lin- 
early upon d and w. In equations (13) through (15), we assume that the particles do not re- 
spond to the turbulent fluctuations, so that the particle velocity distribution function is 
unaffected by the gas fluctuations. From Hinze (1972), this is true if 



(16) 


where A is the integral lengthscale of turbulence, u’ the rms turbulent gas velocity and v the 
kinematic viscosity of the gas. 


Boundary Conditions 


Jenkins (1988) has recently derived boundary conditions for flows of nearly elastic 
but frictional spheres that interact through collisions with a flat frictional wall. Considera- 
tion of the balance of linear momentum, angular momentum, and energy in a single colli- 
sion and some simple averaging provide the rate per unit area at which momentum and en- 
ergy are being supplied to the flow by the wall. Boundary conditions are obtained when 
these are related to the shear stress and energy flux in the flow, evaluated at the wall. When 
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the mean spin of the particles is assumed to be half the vorticity of the mean particle veloc- 
ity and fluctuations in the spin are ignored, these conditions are, in the limit of large fric- 
tion, 

S = - (2 /tt) !/ 2 (N/7w) [v + (d/4) (dv/dr)] ( 17 ) 

and 

q = - (2/it) 1 / 2 N { (l/14w) [v + (d/4) (dv/dr)] 2 - (l-e w +l/7)w } , (18) 

where e w is the coefficient of restitution of the particle-wall collision. Consequently, upon 
adopting the no-slip condition for the gas velocity at the wall and employing the constitutive 
relations (13) and (15), the values of u, v and w at the wall must satisfy 


u = 0 , 

(19) 

dv_ A I 

dr -A d ’ 

(20) 

and 


dw D w 
dr" - B d ’ 

(21) 

where 


A = - 384V2 (1-e) / { 7[20Tt + 96^2 (1-e)] } 

(22) 

and 


B - 175jc {<w )2 ^ ' 14 ^' ew + 7^ ‘ 

(23) 

The second set of boundary conditions is provided by symmetry. 

At the centerline. 

^8* 

ii 

ii 

II 

o 

(24) 

RESULTS AND DISCUSSION 


Equations (1), (2), and (4) are solved numerically using a simple iterative finite dif- 
ference algorithm. To verify the predictions of the model, gas and particle velocity profiles 
are now compared with detailed measurements in a vertical pipe. Using a laser-doppler- 
anemometer, Tsuji, £Lal (1984) measured gas and particle velocity profiles in a vertical 
pipe of 30.5 mm ID. In these experiments, polystyrene particles (p s = 1020 kg/m 3 , d = 
500|im) were suspended in air with ratios of parucle-to-gas mass flow rates (loading) as 
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high as 3.6. Because the integral lengthscale A is of the order of the pipe radius, A/d = 30 
and [(u'A/v) (Ps/p)] 1 ^ ~ 600. Therefore, Hinze's criterion (16) is clearly satisfied. Under 
these conditions, the particles cannot follow the gas turbulence. In the calculation we as- 
sume coefficients of restitution of 0.9 and 0.8 for particle-particle and particle-wall colli- 
sions, respectively. We input the pressure gradient in the gas and the particle concentration 
at the wall and adjust these until we agree with the measured gas velocity at the centerline 
and the measured loading in the experiments. 



Fig- 1- Calculated gas and particle velocity profiles compared with the data from Tsuji, et 
aL (1984). The calculated velocities are normalized with the calculated gas velocity 
at the centerline (7.9 m/s), while the measured velocities (triangles: gas, squares: 
particles) are normalized with the measured gas velocity at the centerline (8.1 m/s). 
The ratio of gas-to-particle mass flow rates is 3.6. 

Figure 1 shows good agreement between the measured mean velocities and the pre- 
dictions of the model. At the centerline, the gas velocity is reproduced to within 5%, and 
the particle slip velocity to within 10%. However, the gas velocity profile near the wall dif- 
fers from the model predictions. This difference is not altogether surprising, considering 
the simplicity of the gas equations used. Nevertheless, the overall agreement in the gas 
phase is good, because the term that dominates the gas momentum equation is not the 
Reynolds stress, but the difference between the pressure gradient and the particle drag. 

The predicted particle velocity at the wall is positive, although not as high as the 
experimental value. Comparable observations were made by Lee and Durst (1982) in a 
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similar situation. At first sight, these observation might be surprising; because the gas ve- 
locity is zero at the wall, one might expect the particles to fall. In fact, there is a region near 
the wall where the particles can acquire a velocity higher than that of the gas. This effect 
results from the shear stress in the particle phase. Particles further from the wall are lifted 
by the gas and, through the particle shear stress, they lift the particles closer to the wall. 
Clearly, the details of the momentum exchange in the particle phase are essential for an ac- 
curate description of the flow. 



Fig. 2 . Particle fluctuation velocity normalized with the calculated gas centerline velocity. 

Figure 2 shows the radial distribution of kinetic energy. In this case, w increases 
away from the wall. There is a slight decrease in energy in the immediate neighborhood of 
the wall followed by an increase into the interior. In more dense flows the supply or dissi- 
pation of energy by the wall is expected to be more important. The particle volume fraction 
is calculated using equation (14), and it decreases away from the wall (Figure 3). Unfortu- 
nately, these predictions cannot be compared with the experiments of Tsuji, et_al. (1984), 
who did not measure the particle velocity fluctuations or the concentration of the particle 
phase. Nevertheless, the predicted concentration profile agrees qualitatively with other ex- 
periments in pneumatic transport (e.g., Boothroyd, 1971 and Kramer & Depew, 1972). 
Several explanations have been proposed for these larger particle concentrations near the 
wall. Boothroyd (1971) attributes this to electrostatic forces, while Berker & Tulig (1986) 
invoke non-gradient turbulent diffusion. In the regime dominated by particle collisions, be- 
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cause of the constant particle normal stress and the constitutive relation (14), this trend is a 
direct result of the profile of fluctuation energy. 



Fig 3. Profile of particle volume fraction. 

Finally, we propose that the enhancement of the gas turbulence observed by Tsuji, 
etah (1984) in the presence of large particles is due to fluctuations in the gas induced by 
the random motion of the particles. Turbulent enhancement of the order of our calculated w 
were measured for the 500pm particles by Tsuji, sUl. (1984), and the strength of the tur- 
bulence was observed to increase with particle diameter. 
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Introduction 

The challenge of the study of gas-solid turbulent flows arises because of the complexity 
of physical interactions between the two phases. Research needs are also not universal for 
different situations. Depending on the particle loading ratios, particle sizes and particle- 
particle collision frequencies, the flow of gas-particle mixture law can be classified as dense 
flows or dilute flows. The characteristic dimensions of the distribution of particles in turbu- 
lent flows may determine whether the two-phase mixture can be regarded as a continuum 
or not [1,2]. In a dilute suspension flow in which particle motion is controlled by the aero- 
dynamic forces on the particle, Crowe [3] has suggested a criterion for treating the particle 
cloud as a continuum. In this case, the Stokes number (St) which is defined St = u r f*/A c 
, where v r is the slip velocity between two phases, is the particle relaxation time and 
A c the distance traveled between collisions, should be less than 0.1 and depending on the 
magnitude of flow Reynolds number, boundary conditions for particulate phase have to be 
modified. 

In this paper, scaling factors determining various aspects of particle-fluid interactions 
and the development of physical models to predict gas-solid turbulent suspension flow fields 
will be discussed based on two-fluid, continua formulation. 

Scaling Rules 

The motions of particles in a turbulent flow field are determined by relative density, 
particle size, inertia, free fall velocity, as well as the correlation between particles and 
underlying flow turbulence. On the other hand, the particulate phase may influence the 
turbulence energy spectrum of the gas phase in wave number ranges corresponding to the 
size of spacing of dispersed- phase dimensions [4]. To investigate the various modes of in- 
teraction, the relaxation time t* of particles has to be compared with various characteristic- 
times of the underlying flow field, t*, in its simplest definition = is a measure of 

how quickly a particle of density p a and diameter d p can respond to changes in the ambient 
fluid velocity, v is the kinematic viscosity of the fluid. In the continuum mixture theory, 
f, can be redefined based on a particle Reynolds number weighted by the concentration of 
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solids [5]. If t* is small compared with a time scale corresponding to some particular flow 
structure, then the particle will follow the motion of that structure; if not, the particle will 
tend to be uncoupled from these motion. 

Hinze [4] pointed out that the dynamic behavior of a discrete particle will not be 
determined by the eddies of size much smaller than that of the particles. The effect of 
these smaller eddies will tend to average out over the particle surface. In the case of 
turbulent flows, this requires that the particle size be smaller than the important dynamic 
scales of turbulence before the time scale ratio becomes a useful measure of the particle- 
fluid interaction. 

Let us consider the dynamically smallest eddies in the turbulent flow. These eddies 

can be characterized by the Kolmogorov microscale r/, and the characteristic time scale r 
1/2 

is of order (u/e) . The ratio of the particle relaxation time to this time then becomes 

d 2 

^ ~ f° r ^he typical suspension problem we are interested in, i.e., p 9 j p > 0( 10 2 ), 

the particles have to be at least one order of magnitude smaller than the Kolmogorov 
length scale in order to be subjected to the motion of the smallest eddies. 

Direct interaction between particles which results from particle-particle collisions can 
be estimated from the ratio of t * and the time scale between particle collisions / c which 
is given [4] ~ 0{ — * 2 — - ) for particles of uniform size d p . v r is the relative velocity between 
particles and n is the particle number density. For the case t*/t c <C 1, the particle has time 
to respond to the local velocity field before the next collision so its motion is dominated 
by the supporting flow forces and the collision which leads to direct interaction between 
particles can be neglected. Then a solid particle is subjected to a variety of time-varying 
forces by the ambient fluid flow. For particles with p s /p > 10 2 , the governing forces due 
to inertia effect ( drag ) and crossing trajectory effect have been singled out [6,7,8]. To 
describe the behavior of particles in a turbulent flow, the simplified Basset, Boussinesq 
and Oseen (BBO) equation has to be solved. This equation in principle cannot be solved 
unless the relation of the Lagrangian and Eulerian correlations of the random fluid field is 
known rigorously. However, the particle trajectory can be determined similar to a random 
walk computation [9] in which a dispersed-phase element is assumed to interact with an 
typical turbulent eddy as long as the relative displacement of the element with respect to 
the eddy does not exceed the characteristic eddy size, / e , and the time of interaction does 
not exceed the characteristic eddy time, t e . The selections of / e and t e are cleanly arbitrary 
since turbulent flows composed a spectrum of length scales and time scales. 

To gain some insight of the scaling rule for turbulent dispersion, the fundamental 
dispersion results of Snyder and Lumley [6] are used to compare with the stochastic pre- 
dictions. The experiments involved the dispersion of individual particles which were isoki- 
netically injected into a grid-generated turbulent flow. The mean flow are uniform in 
the test region and the detailed turbulent structure were measured downstream of the 
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injection point. Since the grid-generated turbulent field is homogeneous decaying and is 
approximately isotropic with fluctuating intensity decreasing in the direction of the mean 
flow. Typically, grid turbulence is also characterized by self-similar spectral distribution 
in which a local set of characteristic eddy time scale can be identified. Following Gosman 
and Ioannides, the eddy time scale is evaluated from the expression t e = / e (^)“ * and 
l e z= Cp k 2 /e, where k is the turbulent kinetic energy (= an< ^ e * s i s °t r °pi c t ur_ 

bulent kinetic energy dissipation rate and C ^ = 0.09. The ratios of t* over eddy time scales 
for the three types of mono-dispersed particles used in Snyder and Lumley’s experimental 
setup are plotted on Figure 1. Mean-squared radial dispersion of the particles is plotted 
as a function of the residence time in the flow for three types of particles. 
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Figure 1. Time scale ratio of various particles in a uniform grid generated turbulent flow 

The results are obtained by averaging over 2000 particles following a 4th-order Runge- 
Kutta integration of the simplified BBO equation ( only inertia term retained ) using 10 
percent of interaction time as integration time step. The agreement between the stochastic 
model predictions and the measurements for the case of corn pollen ( with material density 
lOOOfcg/m 3 and 87//m diameter ) are much better compared to the other two cases of 
hollow glass particles and copper particles. This is not surprising in viewing the scaling 
rules involved ( see Figure 1 ): the corn pollen particles are most closely associated with the 
turbulent time scales responsible for dispersion. For light particle such as glass particles, 
the turbulent eddies responsible for the dispersion should be of higher frequency ( smaller 
time scales ), probably Kolmogorov time scales. The net results is that the numerical 
model underpredicts the dispersion. On the other hand, copper particles should interact 
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with eddies associated with larger time scales. Using the integral time scale (f c ) results in 
overpredicting the turbulent dispersion by inertia effects. 


’a o 



Figure 2. Predicted and Measured particle dispersion in a uniform grid generated turbulent flow 

The numerical models used for the calculation in this study do not account for the 
crossing trajectory effect due to the particle free fall. For the particles studied, especially 
heavy particles, inertia and crossing trajectories are inseparable at the very downstream 
of the decay. According to the scaling rule, gravity has an influence on the two-phase 
flow when n t = u\ where u* is the terminal velocity (= gt * ) and u* is the characteristic 
velocity scale of turbulence. In the grid turbulence jjy- oc {jj) , thus we expect the 

A 1 

free-fall effects to be insignificant in the decay of the grid turbulence if <C 0{ l ~f L ) 1 ^ 2 . 
Here, Um is the mean longitudinal velocity along x-axis and M is the mesh size. Only 
for the smaller particles with t* < 10 msec this condition can be satisfied. However, the 
above calculations indicate the close relationship between the scale parameters and the 
particle-fluid interactions. 

Two-Phase Turbulence Modeling 

In most turbulent multiphase flows of practical interest there exists a spectrum of 
dispersed phase time and length scales. Despite the abundant use of single time and 
length scale models ( such as the k — e model ), the underlying carrying gas turbulence 
is also dominated by a variety of time scales. A two-phase two-scale turbulence model 
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based on continuum approach has since been developed [10] according to the scaling rule 
described above. The modification of particles on the turbulence, the so called modulation 
effect, has been taken into account for large eddies by mean slip between two phase and 
small eddies caused by the particle slip velocity on the fluctuation level. 



Figure 3. Mean axial velocity profiles across the pipe of gas phase and particulate phase 
The turbulent transport equations are summarized here: 
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Here, the k p and k t are turbulent kinetic energy of the large-scale energetic eddies 
and the small-scale transfer eddies and t p is the energy transfer rate from the large eddies 
to the transfer eddies. This model composed of two set of scales for gas-phase turbulence. 
The model is valid for the situation r p ( = r ( Kolmogorov time scale ) and 

particle loading ratio ( p p j p ) of order of 1. This model has been applied to a gas-solid 
suspension pipe flow by Tsuji et al [11]. In Figure 3, the comparisons are made for two 
particle loadings with d p = 200 pm. The predicted velocity profiles are flatter than the 
experimental data and the relative velocities between two phases decrease with increasing 
particle loadings. The distance of the sign change of the slip velocity shifts toward the 
wall for larger particle loadings. The flattening effect by the particles on the fluid velocity 
distributions can be observed in the Figure 4. Besides this flattening effect, the point of 
maximum gas phase velocity even deviates from the pipe axis as the loading increases. 
Such a concave profile indicate counter-diffusion type of momentum transport and cannot 
be predicted by the current model. It is interesting to note that such profiles were not 
found in other LDV measurements for similar configuration [12,13], this phenomena should 
await further experimental confirmation. 



2r/d 

Figure 4. Effect of particle loadings on gas phase velocity profiles 


The modulation effect of particles on the turbulence is shown in Figure 5 for longi- 
tudinal turbulence intensity profiles for 200^/m particles. Cases with particle size greater 
than 200^m were not computed due to model limitation based on the scaling rules. It is 
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known that fluid turbulence is greatly influenced by the particles and the mode of influ- 
ence differs with particle size. For larger particles (> 300 the turbulence intensities 
of the fluid are increased due to the presence of particles while suppression of the turbu- 
lence properties is observed for smaller particles. It is seen in Figure 5 that for a particle 
loading of 0.9 the turbulence intensities are reduced by 30 % in the core region. However, 
the intensity in the core region increases again as the particle loading increases from 0.9 
to 3.2 and the intensity in the wall region is monotonously damped. This phenomenon 
can be explained by the competitive mechanism between modulation due to mean motion 
and fluctuation motion. The small scale modulation effect always acts as a sink term in 
the k t equation while the large eddy motion modulation effect can be extra production 
or dissipation depending on the signs of (Vi ~ V t ) and the distribution of mean particle 
density profiles. The cross-over of the intensity profiles is closely related to the cross-over 
of mean velocity profiles of gas phase and particulate phase. The relative magnitude of 
the modulation effects is reasonably well represented by the model. 



Figure 5. Modulation Of particles on gas phase turbulence intensities 


Summary 

The modes of particle-fluid is discussed based on the length and time scale ratio, which 
depends on the properties of the particles and the characteristics of the flow turbulence. 
For particle size smaller than or comparable with the Kolmogorov length scale and con- 
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centration low enough for neglecting direct particle-particle interaction, scaling rules can 
be established in various parameter ranges. The various particle-fluid interaction will give 
rise to additional mechanisms which will affect the fluid mechanics of the conveying gas 
phase. These extra mechanisms are incorporated into turbulence modeling method based 
on the scaling rules. A multiple-scale two-phase turbulence model has been developed, 
which gives reasonable predictions for dilute suspension flow. Much works is yet to be 
done to account for the poly-dispersed effects and the extension to dense suspension flows. 
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ABSTRACT 

This paper presents a multiphase turbulence 
closure employing one transport equation, namely 
the turbulence kinetic energy equation. The 
proposed form of this equation is different from 
the earlier formulations in some aspects. The 
power spectrum of the carrier fluid is divided 
into two regions, which interact in different ways 
and at different rates with the suspended 
particles as a function of the particle-eddy size 
ratio and density ratio. The length scale is 
described algebraically. A mass/time averaging 
procedure for the momentum and kinetic energy 
equations is adopted. The resulting turbulence 
correlations are modeled under less restrictive 
assumptions comparative to the previous work. The 
closures for the momentum and kinetic energy 
equations are given. Comparisons of the 

predictions with experimental results on 
liquid-solid jet and gas-solid pipe flow show 
satisfactory agreement. 

NOMENCLATURE 


a amplitude ratio 

b body force 

C -’ • C *2' C *3' < >r c *s- S' 


'V *1 ' 


C = constants 
P 


d 

D 

E(k) 

J. 


particle diameter 

turbulent diffusion coefficient for solid 
phase 

energy spectrum 
flow variable 
interaction term 

flux vector for a variable $ 


k kinetic energy of turbulence 

£ length scale 

Ng Stokes number 

K(r , t ) phase distribution function, Eq (1) 

p pressure 

r posi t ion vec tor 

Re Reynolds number based on the most 

energetic eddy size 

t.AT.At time and averaging time intervals 
corresponding to the turbulence 


( ) - For further correspondence; 

On sabbatical leave at Caltech 104-44, 
CA 91125, until lay 30, 1989. 


V , Av 

X . ,x . 
1 J 

y 

a 

t 


V 

*C 

V 
p 
a 

T 


production and transfer range, 
respectively 

Eulerian and Lagrangian velocities, 
respec t ive ly 

volume and averaging volume, respectively 

x Cartesian coordinates 
n 

distance from the wall 

concentration of volume 

dissipation rate of turbulence kinetic 

energy 

Kolmogorov length scale 
wave number 
dynamic viscosity 
densi ty 

turbulent Schmid t/P rand t 1 number 
shear stress 


Subscripts 

i,j,n denote Cartesian coordinates (= 1.2.3) 
e eddy 

k turbulence kinetic energy 

K flow component K 

KP phase K in the production range 

KT phase K in the transfer range 

£ laminar 

L.S denote liquid and solid, respectively 
P production range 

t turbulent 

T transfer range 

Tot total 

2 average over area 

t dissipation rate of turbulence kinetic 

energy 

Superscripts 


f mass average of f 

f* turbulence fluctuation of f 

{• fluctuation of f at low wave number 

(production range) 

f fluctuation of f at intermediate wave 

number (transfer range) 

P production 

T transfer 

7? { v volume/time and mass/time averages of f 
K 

over K 

<f > intrinsic space average of f over K 

K 
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1 • INTRODUCTION 

Multiphase flows are widely applied In 
engineering processes from chemical, petroleum, 
mining and other industries. Various theoretical 
and experimental techniques for the investigation 
of those flows are available. Some of them are a 
straightforward extension from the single phase 
flow models by introducing some ad hoc 
modifications. Other investigations originate 
from the gas-solid flow (Soo. 1983) or fluidized 
bed models (Wang et al . , 198S). 

Increasing concern for the prediction of 
turbulent multiphase flows have been noticed 
during the last twenty years (Danon et al . (1974), 
Al-Taweel and Landau (1977). Genchev and Karpuzov 
(1980) . Melville and Bray (1979). Crowe and 

Sharma (1978). Michaelides and Farmer (1983). and 
Shuen et al . (1983)). Two equation turbulence 

models have been proposed for dilute particulate 
flows by Elghobashi et al . (1982. 1983, 1984) and 
Crowder et al. (1984). Algebraic and one equation 
turbulence models have been suggested also for 
dense liquid-solid flows (Roco et al . . 1983. 1985. 
1986) in which the particle-particle interactions 
play an important role besides the fluid-fluid and 
fluid-solid interactions. Most of these studies 
as well as other earlier investigations have some 
limitations. In the above mentioned studies the 
response of solids to the turbulent fluctuations 
of the carrier fluid is obtained under 
restrictions similar to those refered by Hinze 
(1975. p. 460). which limit their use. In 
addition to that. empirical constants and 
empirical functions are usually introduced in 
these models. 

The purpose of the present paper is 

i) To propose a specific mass/time averaging 
approach for multiphase turbulent flows. Even if 
the approach is developed for incompressible flow, 
its application for other multiphase flows is 
foreseen: From the liquid-solid interaction forces 
only the drag force Is considered in this paper. 

11) To liq>rove the one-equation turbulence 
■odel reported In [28] by including the modulation 
of turbulence by particles as a function of 
particle size and density. 

iii) To test the proposed model with other 
models and experimental data for various two-phase 
flows, without adopting any adjusting empirical 
coefficients. 

2. MIXED AVERAGING APPROACH 

The continuum transport equations for 

multiphase flows can be obtained by assuming a 
continuum medium with averaged field quatities by 
using either time, local volume, local mass or 
spectral averaging (see Buyevich (1971). Soo 

(1967), Vernier and Delhaye (1968). He tsroni 
(1982)). The averaging for multiphase flow 
systems may be performed in various ways. Mass 
averaging technique was applied by Abou-Arab 

(1985) for turbulent incompressible and 
compressible flows. To express the spatial 
nonuniformities and interactions between the flow 
components Roco and Shook (1985) have developed a 
specific volume/time averaging technique for 

turbulent mul ti component systems, in which the 
size of the averaging volume Av is related to the 


of rh#» r?* S f aie * Since the Eulerian description 
dMrr<nM° W * t J irc convenlent than the Lagrangian 
fwthpmr<° n | t * lere are more comprehens i ve 

schem «s for such formulation, they 
tin*. ° rme V he V0lu,ne/tlme averaging into double 
>g r ng ' For flow Action f and any 

component K in the mixture, at a position r and 

time t 


fK (L. t) 

where 


1 

AT 


c 


•t+AT/2 


■AT/2 


< f K > dt 


( 1 ) 


AT = time averaging interval corresponding to 
the turbulence production range (AT ®) . 
j rt+At/2 

f K > = "77 I *(£.t) K(r.T) R(T-t) dr 

J t-At/2 

(2a) 

* Intrinsic averaging of f over At and 
the flow component K. 

At = Eulerian time scale for the most energetic 

eddies =J R(r-t)dT (corresponding to the 
0 

Taylor’s length scale). Note that At<<AT 

but much larger than particle residence time 
at £ . 

1 if the considered flow component 
K(£.t) = | resides at point r^ and time r . 

0 otherwise 

= phase distribution function 


(2c) 


R(T-t) = V "( 1 ) v ‘*f t-r) 

> 2 (t) 

= auto-correlation coefficient of the 
velocity fluctuations of the most 
energetic eddies (Taylor s scale) 

The time averaging over At corresponds to the 
averaging over local volume Av. The dimension of 
y * 8 lven by the turbulence mixing length, and 
= (length scale of Av)/(mean velocity) (291 
By Integrating over a flow component K its 
aC ** f " lth other flow components become 
boundary of Integration, and the interaction terms 
,. r , derlved In a straighforward manner in the 
differential formulation. According to (1) anv 
Instantaneous value differs from the mean value by 
a turbulence fluctuation with two components f. 
and f": K 


f K - f K + ft = f K + f K + f K 
where 


(3) 


f“ - t 

K ~ f K 


< f K > 


- spatial nonuniformities within Av (or 
temporal nonunlformi ties within At) 


f K = < f K > " ** 


<■*) 


(5) 
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s temporal nonuni formi ties of < > within AT. 

Since the averaging domain (Av or At) has the 
dimensions of the mixing length or Eulerian time 
scale, respectively, the turbulence fluctuation fj^ 

corresponds to the turbulence transfer range. The 
temporal nonuniformities f^, reflect the turbulence 

fluctuations in the production range. 

By averaging with formula (1) the point 
ins tantaneous conservation equations, one obtains 
the double time averaged equations. The 

formulation is equivalent to the volume/time 
averaging. The momentum and kinetic energy 

equations are given in Appendix A. 

The local mass average of f over a flow 

component K is denoted f^. It is obtained by 

applying (1). in which the phase distribution 
function K( - t,t) is weighted by the specific mass 

v 

In this paper we model the phase interaction 
by using spectral analysis and suggest a closure 
of the mass averaged equations for linear momentum 
and kinetic energy. The averaged equations are 
initially written with all the terms, and then 
simplified formulations for various flow 
conditions are suggested. 

3. ENERGY SPECTRUM AND SOLIDS - EDDY INTERACTION 

It Is well accepted that turbulence is 
characterized by fluctuating motions defined by 
an energy spectrum (Tennekes and Lumely (1972)). 
Single time scale models, which are normally used 
for the prediction of turbulent flows, seems 
simplistic because different turbulent 

interactions are associated with different parts 
of the energy spectrum (Hanjalic’ et al . (1979)). 
A typical energy spectrum can be divided into 
three regions, The first region is the production 
region of large eddies and low wave number. The 


third region is the dissipation region with small 
eddies and high wave number. In which the total 
kinetic energy produced at the lower wave number 

is dissipated. The intermediate range of wave 

numbers represents the Taylor's transfer range. 
The total kinetic energy k of turbulence may be 
divided into production range (k p ) and transfer 

range 

(Itj.) because there is negligible 

kinetic 

energy 

in the dissipation range: 


where 

k — kp + Jcj. 

(6) 

k P 

1 ‘2 
= 2 u i ■ 

(7a) 

4 

1 ~"2 
= 2 U i 

(7b) 

V 

, = fluctuating velocities in the 



production and transfer range, 
respectively . 




Wove Number, /c(m'' ) 

Figure 1. Schematic showing the relative particle 
size to different eddy sizes in the 
energy spectrum. 


This partitioning of the energy spectrum was shown 
to be important for swirling flows (Chen (1986)). 
and hetrogeneous mixture flows such as two-phase 
Jet (Al-Taweel and Landau (1977)). 

By using spectral analysis in conjunction 
with mass/time averaging some additional turbulent 
correlations will result in the mixture flow 
equations comparative to homogenous flows. These 
correlations can be classified into five 
categories : 

1) Eddy-eddy interaction 
11) Eddy-mean flow interaction 

iii) Eddy-particle interaction 

iv) Particle-mean flow interaction 

v) Particle-particle interaction (for dense 
suspension flow). 

These correlations have to be modeled. Since the 
suspended particles may be of different sizes and 
different materials, their response to the carrier 
fluid fluctuations will vary as a function of the 
mean and fluctuating properties of the flow. The 
present work will consider a two-way interaction 
mechanism between solid particles and fluid 
vortlcies in dilute suspensions. This interaction 
mechanism depends on the ratio between the 
particle size dg and the turbulent vortex (eddy) 

size 1^. These length scales are compared with 

the turbulence dissipation micro-scale (q) . 

To determine the particle-eddy Inter action 
the energy spectrum for multiphase flow system is 
divided Into three typical zones (Figure 1): 


1. "Large vortex zone" (#1). where the 
turbulence energy is extracted from the mean 
flow by low frequency eddies Here, the eddy 
length scale € is larger than the particle 

size d s : 
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2 . 


‘el > d S > ” 


(8a) 


"Medium vortex zone" (82). where the solid 
particles are about the same size with the 
vortex size . 1 . e . 

*e2 ~ d S > 11 (8b) 

"Small vortex zone" (#3). which would 
correspond to the Kolmogorov's length scale, 
i . e . 

d S > ! e3 * *> (8c) 


+ U Ki Vkj + “KjViii + “Mj + VKi u Kj 

+ “K^i^i + ^K^h'j + VKj u Ki + a K u Ki u Kj ) 

= P K^“K b Ki + “K^Ki + “?Ki + Sjfti + ^Kil 
Body force 

d r -~ — , 

ax. Lq k p k + “kPk + “k p k + “k p k + 

Pressure effect 


In zone #1 the solid particles generally follow 
the motion within a vortex, and have an energy 7 
dissipation effect. The particle response to the 
turbulent fluctuations (turbulence modulation) is 
fully determined (see Hinze, 1975). In the snail 
vortex zone #3 the solid particles can not 
significantly affect the turbul ence 
microstructure. For the intermediate zone #2 a 
linear variation of the particle response is 
considered. This petitioning allows for the 
par t 1 c 1 es~nonun 1 f orm size eddies interaction to be 
efficiently modeled. 

The present closure formulation originates 
from the idea of subgrid scale modeling. If this 
Idea is to be accepted, any flow quantity u. v, a. 

^ etc * may be separated into three parts 

according to (3). where f R and f" define the 

fluctuations in the production and transfer ranges 
of the energy spectrum, respectively. By starting 
from the particle equation of motion in Its 
general form, the relation between the particle 
motion and different fluid eddies can be 
determined, and from here the fluctuation 



Frictional effect 


+ “iPibji + Viaid ♦ n Kt > k O) 

Riase interaction 


where K is a flow component. 

in the i-th direction. 

projection in the i-th 
interaction vector (I ) 

-K - K 


is the body force 

a ™ 1 (I Ki^-K ls the 
direction of the 


Equation (9) contains correlations which are 
related to the production and transfer wave number 
ranges of the turbulence spectrum. It contains 
also mixed correlat ions . 

^ ^ The Kinetic Energy Equation 


components f^ and f£ (see section 6). 

4. COMPOSED AVERAGED EQUATIONS 

4.1 Mean Flow Governing Equations 

The mass/time averaged momentum equation 
(Appendix A. Eq. A2) with f* = + f" yields: 

P K 3? ( “K u Ki + “k“k + “k“k + a u Ki + “"“Ki) 

Time rate change of the mean fieri? coru>ectlon 

+P K £- t“K\i“Kj3 
J J 

Mean flo® connection 


+ P K 3^ [ “K u Ki u Kj + °K u Ki U Kj^ 
Inertial effect 


+P K ^ a K u Ki u Kj +a K u Ki u Kj^ + u KjK u Ki +a K u Ki^ 
Col l i sional/Inert ial Effects 


+ U Ki t“K u Kj + Vi<j] + a K^ u Ki u Kj + U Ki U Kj^ 


Similar to the mean flow governing equation 
the mass/time averaging form of the turbulent 
kinetic energy equation can be derived for a flow 
component K. The exact form of this equation for 
steady state turbulent flow is given in Appendix A 
(Eq. A5) . It contains more than one hundred 
correlations which are related to the eddies in 
the production and transfer ranges as well as some 
mixed correlations. However, only some of these 
are predominant in a given flow situation as a 
function of relative particle-vortex size and 
density ratio. 

* ' ^ S ome Modeling P rinciples and Assumptions 

By analysing the derived mean momentum and 
kinetic energy conservation equations, it can be 
easily recognized that some modeling assumptions 
must be made, based on the physical interpretation 
and the nature of each term. Previous 
experimental and theoretical findings can help in 
modeling collectively similar terms with minimum 
number of empirical constants. Secondly, carrying 
out an order of magnitude analysis for different 
correlations which appear in the governing 
equations some terms may be neglected. Thirdly, 
the micromechanics which control the ability of 
the flow variables to correlate with each other 
and the factors affecting the magnitude of these 
correlations should be considered. With the 
previous remarks In mind, one can assume for the 
sake of simplicity that: 
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1) The correlations between the large eddies 
from the production range and the snail 
eddies from the transfer range (mixed 

correlations e.g. u’u") can be neglected as 
they originate differently and they are 
related to different ranges of the power 
spectrum. Similar assumptions are accepted 
in the classical single fluid turbulence 
theory. 


where 

Ap = C fiP p L k LP V 

At = C jjT p L k LT 1-j- . and (16) 

At = Ar + At = C M p L k L 1/2 A (17) 


The void fraction fluctuations occur mainly 
at low frequencies i.e. 


a = a + a * + a" 


( 10 ) 


with a* >> a” 


(in 


This is a simplifying assumption which is 
acceptable for such complicated problems. If 
the particles are of small diameter their 
concentration is relatively uniform 
distributed in the Taylor length scale. High 
void fractions are mostly associated with 
large size and high density particulate 
flows. These large size particles are mainly 
fluctuating at low frequencies due to its 
high inertia, hence they in turn correlate 
weakly at high frequencies. 


3 ) 


The correlations of higher order than three. 


for instance a’u’dp'/dXj, 


. . dp_ 
Vj 


etc., are neglected. These are at least an 
order magnitude smaller than those of the 
third order (see Hanjalic’ and Launder 


(1972)). 


4) Pressure diffusion contribution to the total 
turbulent diffusion in the kinetic energy 
equation will be neglected because of its 
relatively small magnitude (Hanjalic’ and 
Launder (1972)). 


5) The Boussinesq gradient type approximation is 
adopted for modeling of different fluxes and 
triple correlations, with assumptions similar 
to Elghobashi and Abou-Arab (19S3) and Roco 
and Mahadevan ( 1986) 


6) The following constitutive relations are 
employed for the shear stress of carrier 
fluid: 


chjj du . 

p Vj = Ap ( 5^ + 


3 k LP 5 ij 3 AP 6 lj u n.n 

( 12 ) 


- pu lUj 


*1 


cKi 


cbc , 


' 3 k LT 6 ij 


^LT^i j U n . n 
(13) 


Similar relations can be written for the viscosity 
of the dispersed phase (see Roco and 

Balakrishnan (1985)). However, in the present 
work we choose to define the eddy viscosity of 
solids as follows: 


Dp = d, /a 
St Lt 


where 
and 


O - o / o 

aS S 


°aS = u Lt^°S 


a S = "st^S 


(IS) 

(19) 

( 20 ) 
( 21 ) 


Appropriate expressions for o ^ are cited in many 

articles such as Peskin (1971), Picart et al . 
(1986), and Hetsroni (1982). is a Schmidt 

number and its value 5s about 1.5 (Abou-Ellail and 
Abou-Arab (1985)). 


7) In the present approach for dilute 

particulate flows the turbulence kinetic 
energy equation Is written only for the 
carrier flow. The solid phase turbulence 
kinetic energy and turbulence correlations 
are evaluated from the available solution of 
the linearized equation of motion of a solid 
particle after its transformation from the 
real time to the frequency domain. 

8) Terms which are of similar nature i.e. 

convection, diffusion, d i ss ipat ion . etc. can 
be modeled collectively. The length-, 

velocity- and time-scale which are 

appropriate for the description of their rate 
of change must be identified from the 
physical interpretation of these terms. 

9) The response function which shows the ability 
of solid particles to follow the eddies is 
obtained from the equation of motion of 
particles for different local dimensionless 
parameter , dg/1^ and Pg/p^ . 

10) To establish the degree of generality of the 
proposed model validation tests were carried 
out for air-laden and water-laden jet. and 
aii — solid pipe flow. 


5. CLOSURE FOR THE MEAN FLOW EQUATIONS 


while the total shear stress -p u.u^ 03X1 ^ gl yen 
by 


P 


t t 
u . u . 
i J 


p u ! u '. 
i J 


p u'.'u'.' . 
1 J 


(H) 


With the modeling assumptions given in the 
previous section, the steady state mean flow 
momentum equations for any flow component, reads 


d - - 


p k ET S (‘VkiV + dT s p k ( ‘YvAj + °K u Ki u kj ) 
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+ 4 “ Ki (VKj 4 

4 "Kj (^Ki 4 %)] + (Vtj, - ^iaji)} 

= p A\i + °K b Ki + “K b K i ^ ~ 4 ^ + “k p K ) 


4 (P, 


- *K 


K &x 


i 


Pk ax. 



The diffusion fluxes namely a^u^ . , . . . etc. 

are modeled using Boussinesq approximation as: 


-r—r _ U KP ^“K 

VK1 ~ P &r 

°oK 1 


(26) 


while 

- U KT ^ 

V'Ki T dx~ 
°aK 1 


(27) 


+ 3Xj ( a K T iaji) + ^Ki^-K f 22 ) 

In Eq. (22) there are 16 correlations, half 
of them in the production "large eddy" range of 
the spectrum. The terms are modeled following the 
criteria- i) Physically correct behavior, ii) 
Minimum number of empirical functions and 
constants. and ill) Comparisons against 
experimental data over a wide range of conditions 
are required to check the validity of the model. 


Similar expression can be written for aLu’ and 


^K^KJ ’ ' * ' etc - However, according to our 

assumptions, especially those concerning the void 
fraction fluctuation at high frequency a", the 
fluxes in Eqs . (26) and (27) can be modeled 
col lec t ively as : 


XT . ^ Kt^K 

Ki ~ o v dx 
aK i 


(28) 


The turbulent stresses caused by the large 
and small size energetic eddies. -p^ujuj and 

~ P L U i U j * are defined by Eqs. (12) and (13). The 

correlation between u" and uj* is weak when i / j. 

This finding will be explained in the next 
sec t Ion . 

The colli sional effect correlations are 
modeled after Launder (1976): 

\p 


WKj = “ C *5 ^ {u Ki u Kl d^ u Kj a K 


+ U Kj U Kl “Ki 0 * 5 


(23) 


where is a constant of a value of about 0.1 


Similarly 




~ c ^5 (u Ki\l axj U Kj°K 


The solution of the two transport equations 
(for kp and kp) would require also a description 

for the length scales. This point will be 

explained in more details in the next section. 

The fourth group of terms with the void 
fraction-shear stress turbulent correlations 
contains the laminar viscosity as a multiplier, 
and will be neglected due to its smaller order of 
magni tude . 

The pressure effect contribution to the mean 
flow equation consists of two groups. Each 
contains three terms. The first of these terms is 
the mean pressure-void fraction. The second and 
the third terms are the pressure-void fraction 
correlations which can be modeled after Elghobashi 
and Abou-Arab (1983) by 

" <=** + W = h 4 *2 (29) 

where 


. _ r ,1/2 t t 

h C *3 P X Sc ‘ “fcA 
and 


U Kj u Kl 


a 

dx , 


U K i ct K ^ 


(24) 


C f4 p K k 


1/2 


t t 

“X/A 


where must be optimized by comparison with 

experimental data The fluxes in (23) and (24) 
can also be collectively modeled as: 


a K u Ki u Kj 4 VKiSij = 

^ K.tt a tt 

♦5 ' U K i °X 1 ' dXj \j“K 

+ “Kj^l ' dx]" u Ki a K ) f 25 ) 


The values of the constants C and C are about 

T ^ 

uni ty . 

The second correlation in the second group of 
terms can be also modeled following Launder (1976). 
The final form is 


**X dx. 4 Pjc dXj " P K t C fl‘ U Ki a K 
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t t 

u.,.u, 


„ , Ki K1 2 £ x t 

C *2 ( " 3 1 * U K1°K 3 


The first group of terms in the k-equation 
(Group #1) is the convee t ion of the total specific 
kinetic energy, where 


du 


du, 


+ P K (° 8 u ki“k ' ftT 1 " 0 2 u Kl“K ' dx K1 } (30) 


dx . 
i 


The values of the constants and C *2 are 4.3 and 

-3.2. respectively. The modeling of these 
correlations suffers from the embodied assumptions 
concerning the velocity and length scale 
description. It would require a large number of 
transport equations to model accurately each of the 
above correlations. 


The interaction term (I^)_k ^ or ^ = L 

(liquid) and -K = S (solid) is modeled for dilute 
suspensions with particle Reynolds numbers less 
than unity- 

Th~h- - < 1 S V d S> - “Si> “s 


A p 
♦ (— 

-LT , *L # 
+ — > dT * < 

V L? 

v ll . 
+ T > 

da 

*7) 

? 

a oS 

°aS 1 

°°aS 

°°aS 

where 

the first term 

is 

the drag 

interaction for 


particles in the Stokes range. The second and 
third terms are the turbulent fluxes due to the 
relative motion between the particles and fluid. 
The gradient transport model with the exchange 
coefficients and corresponding to the 

production and transfer ranges is adopted for 

these fluxes ( . o£ . a $ ujj. and 


ax . 

j 


ak 


LP 


ak 


LT 


dx . 
J 


dx t 


(33) 


Diffusion transport of k is composed of two 
main parts. The first part (Group #2) is the 
velocity diffusion and it contains the 3rd order 
velocity correlations, while the second part (Group 
#3) is the pressure diffusion with the 
pressure-velocity correlations. The modeling of 
the velocity diffusion part is obtained as follows 


- — ; — - U LP ^LP 

P L “L u Lj u Li u Li = _ p l“l ~ dT~ 

a k J 

and 

- — — ^ — — - "lt 

P L °L U LJ U L1 U L1 = - P L“L — ET- 

a k J 


(3-1) 


(35) 


P T 

where 0, and o, are the turbulent Prand 1 1/Schmi d t 
k k 

numbers for the kinetic energy in the production 
and transfer ranges. To reduce the number of 
empirical constants and the number of governing 
equations the above two correlations are modeled 
collectively as follows 


- p L u LJ u Li u Li 


- Lt 


ak, 


' p l“l 


0 . dx . 
k J 


(36) 


where is of order of unity. 

The presure diffusion term is negligible 
relative to the velocity diffusion (Hanjalic* and 
Launder (1972)). 


If only single velocity scale is chosen for 
the whole energy spectrum, k^, there will be only 


one 

momentum 

exchange 

^LP 

and i> LT 


U Lt 

- ULP 4 

■lt 

P 

T 


00 0 
aS 

00 o 

aS 


instead of 


(32) 


6. CLOSURE FOR THE KINETIC ENERGY EQUATION 


Mixed and higher order correlations (Croups #3 
and #5} can be neglected according to the modeling 
assumptions stated and discussed in section 4.3 of 
the present paper. 

The production terms are divided into two 
groups. The first group (Group #6) is conmon for 
single and multiphase incompressible flows. 

^Li - —r. — :r- 

P L“L u Li u LJ • ET~ * p L“l u Li u Lj • 
be modeled collectively as follows' 


du. . 

Li 

dx . 

J 


and can 


In the present work, the turbulence kinetic 
energy for the liquid phase **k^" (turbulence 

velocity scale) is obtained from an exact transport 
equation (Eq. A5 in Appendix A), and the length 
scale “1" is described algebraically. The kinetic 
energy equation A5 contains a large number of 
turbulence correlations. In order to obtain an 
engineering turbulence model, it is sufficient to 
consider the principles and the assumptions given 
in Section 4.3. By engineering turbulence model, 
it is meant, a physically correct model with 
minimum number of empirical coefficients. 


— t t 

- 3 u 



Ai 2 

( — Li 

[ dx. } 

(37) 

p. a, u. . u. . • _ 

L L Li Lj dx . 

J 

V\t 

The 

physics 

of 

turbu 1 ence 

and the 


consideration of the spectral energy transfer 
assume that the production is only due to the 
interaction between the mean flow and the large 

eddy. Since the uj^u”^ correlation is for medium 

size eddies which have almost no direct interaction 
with the mean flow, it results that the 
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contribution of these small eddies to the 

turbulence production via the mean field is smaller 
than that of the large eddies. This 

means also that correlation Is weak If 1 / 

J, i.e. It is only of significant value if the 

2 

turbulent normal stress components (uj^ . 1 = 1, 2, 

3) are considered. According to Hanjalic’ et al. 
(1979) multiple scale model the turbulent viscosity 
is defined as follows: 


du , 


Terms like ~ and are modeled in 

a similar manner to that of the above terms i.e. 


du 


' (a L u Lj + Vlj) 



dx . 

l 


U Lt 

°aS 


ax 4 ax. 


(45) 


The £ * tra — Production terms (Group 87) can be 
written in the following form 


= C M k L ( k LP /t LP^ 

(38) 

l L = k LP + k LT 

(39a) 

= 0.09 

(39b) 


This equation can be rewritten as 


v dcL dp du. . 

Extra production = -=^- — -E ( — - + u ki 

CT aS **1 «*! L J 3Xj 


u , r 3a, _ 5u 

^KT^isr 1 ) M6) 

V s J J 


^t /p L = C u 


■ k LP + k LT k LP 
fc LP 6 LP 


= Vi?\? * VlT^LT " ^P + h.T> /p L W 

where C _ and C _ are two additional constants and 

pP pT 

l^p and l^p are also two additional length scales 

for the large and medium size eddies. The length 
scale can be related with the following relations 
(from Eqs. (3 S) and (40)) 

*LT = C 1 J LP ^LT^LP^ 

Since the ratio k^yk^p is of the order of unity 

and 1 > 1 the constant should be smaller 

than unity. Thus If the multiple time scale model 
is not recommended (due to its large number of 
additional constants) an alternative approach is to 
consider a multiple velocity scale model. In this 
model only two differential equations for k^p and 

have to be solved. The length scales can be 

obtained by using the previous relation Eq . (41) 

and any expression for the length scale of the 
large eddies l^p, for example that used by Roco and 

Shook (1983) for cylindrical pipes. 



dissipation is mainly confined to snail scale 

eddies and to simplify the mathematical form of the 
nxiltipie time scale turbulence models, Hanjalic’ et 
al. (1979) have, assumed that there is an 

equilibrium spectrum energy transfer between the 

dissipation and transfer region i.e, e - e 

' total ~ L' 

where t)L Is the dissipation rate in the single 
scale scheme. Thus 


"Wlp = (47) 


" P L°L t LT (4S) 


^Ll ^Ll 

WLarsr 

J J 

and 


- ^Li ^Li 
J J 


where 


£ LP " Sp k Lp 5/1 LP ^ 

£ LT = Srr k LT /1 LT 


The modeling of the additional production terms 
(Group 87) is achieved as follows: 


— r — r— V 

^^1 dx . P dx, dx. 

1 ° aS 1 1 


(42) 


°L U Li dx T dx dx 
1 CT oS 1 1 

and both collectively as 


(43) 


u L t **1 

i aS 1 1 


(44) 


In equation (47) a spectral 
production and transfer 
(Hanjalic’ et al . (1979)). 
additional relation between 


cascading between the 
eddies is considered 
This equation gives an 
the spectrum scales. 


The correlations between the fluctuating 
velocity component and the fluctuating friction 
orces (Interaction terms in Groups #8 and #10) are 
due to fluid-fluid and fluid-solid drag force in 
ilute flows. The friction interaction terms due 
to molecular collision (fluid-fluid interaction. 
Group #8) are given above by Eq . (47) and (4S) 

The form of the correlation between the fluctuating 
drag force and the velocity fluctuations depends on 
the expression adopted for the drag force The 
viscous drag correlation (VDC) in Group #10 for 
Stokes flow over particles in dilute suspensions 
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with the following dimensionless parameters: 


18m, 

VDC = 2~ (U L1 " U Si )-A l 

d S 


18 h - 

- —T (a s' A 2 + A 3 5 


s 

where ‘ 


*1 = a s u i 4 a s u i 

A 2 = u Li (u Li" U Si 1 4 U Ll ( “Lr U Si ) 


(49a) 


(49b) 

(49c) 


s = P S /P L 


N S = 


J U L Kd| ^ 


where 


xRe 

*T 


tCp = Is the wavenumber of the most energetic 
eddies , 

Re = Reynolds number based on ly. and 
Ng = Stokes number. 


A 3 = a S u Li (u Li- U Si ) 4 VLi (u Lr U Si) (49d) 

The first correlation group ^ (Eq. 49a) can be 

approximated using the gradient type assumption. 
The second correlation in this expression (49b) is 
that due to the relative slip fluctuating motion 

u U< u U - U Si> and U Ll (u Li - U Si } - 111656 6311 1,6 
modeled using similar approach to that of 
Elghobashi and Abou-Arab (1983) but with some 
modifications which allow for different 
particle-eddy interaction according to their 
relative size. These modifications are based on 
the spectral analysis carried out by Dingguo (1987) 
for the response of the particles to the turbulent 
fluctuations of the carrier fluid. For very large 
eddies x << Xg. 


For small eddies with wave numbers x >> x g . 
the particle response can be described by 


*S 3 P L 

u s = < V sh - ( V L>k Hf> ^ 

For the intermediate size eddies x 
can use (50) with 


(53) 

- *s one 


a = 

P = 


1_ 

s 


Bind 


-1 

tan 


(s) 


(54a) 

(54b) 


If the fluctuating slip velocity is defined 


as 


w . 

l 


Si 


Li 


(55) 


Ug = ( V $) K * (\) K a * cx P[“i(* u L t-0)] (50) 

where x is the wave number; Xg is the 

Basse t-Boussinesq-Ossen wave number defined as 
l/d^\ (Vg) and (v^)^ are t ^ ie so ^ icl liquid 

velocity components with the wave number x; a is 
the amplitude ratio of oscillations 

» = [(!♦ q 2 ) 2 ♦ q 2 ] 0 ' 5 ( 51a > 

and p is the phase angle of oscillation 

P = tg 1 [q 2 / ( 1 + Rj)] (51b) 

The expressions of and are 


~2 , . 2 

then the ratio of the mean square w. and u^ 


becomes 


4 “L?> 7 U Li 


wi th 



1 , 2 

(1 ♦ 

• 2 . 

u Li U Sl 

= 2 U Li 

U Si ' 


1 "T 2 f . 
= 2 U Li (1 

4 r SL 

‘ r RlJ 

The values 

of r sL and 

can be , 


• 2 ' ' 2 - w: 2 /u,h 


(57) 


the preceding solution (Eq, 50)' 


9N S i_ s 

(1 + — — )(”o%) 
l2( s+0.5) ° 


SI 


(s+0.5) 


2 2 

2 <2*S + ^ 4 (1 4 


9N 


(52a) 


S ,2 


■f2( s+0.5) 


N, 


9 < l (2Hg + /) 

(s+0.5) 2 S 2 

q-. = 1 FIT 9F 


(52b) 


^-2 i*i 4 J 2 4(1 + -^- )2 


(s+O.5) 


l2( s+0.5) 


r gL =[J a^f^d* + j” {-|) 6 s' 2 E L (K)d K ]/J^E L ( K )dK 


r SL = tan ’ 1( V K T ) 

where 

ISu 

h = — (Kg/^) 2 ^. A = 


(58) 

(59) 


k-j.u l (s+0. 5) 


(s + 0.5)d 
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Figure 2. The effect of eddy-particle size ratio 
(* 5 ^* 7 ) on t ^ e particle response to the 
-“2 — 2 

eddy fluctuations ) 


and 


r RL = (» 2 -l) E l («) d* 

+ f (( T )3 s " Ll ) 2 / I E(k) d* (60) 

K s 0 

where aj = (f 2 + f 2 )^ and E^(k) is the liquid 

energy spectrum function for which any arbitrary 
form can be adopted, for instance 


.2 

2 U Li 




“T H 




(61) 


Figure 

eddy-part icle 


2 illustrates the effect of 
size ratio expressed as Kg/*-]- on the 


ratio 

response 


2 2 

u' /u’ which 
Si Li 

to the eddying 


indicates 

motion 


the particle 
1 1 can 


be noticed that high values of (l.e. snail 

particle or large eddy) the particles follows quite 
well the eddy motion. 


Substituting the expressions for r and T D1 

oL KL 

(Eqs. (58) or (59), and Eq. (60)) into Eq. (57). 
the correlation 4„ can be obtained 

A 2 = 2 U Li " f SL + f RL } (62) 


The above analysis applies equally well to the 
large eddies as to small eddies. The energy 
spectrum function E^(k) for a two-phase flow is 

given by Al-Taweel and Landau (1977). However, 

since its form is not essential (Dingguo (1967)) it 
is sufficient to adopt any simple form as that 
given above by Eq. (61). It is clear from the 
above analysis that the particle response to the 
carrier fluid fluctuations is a function of the 
density ratio pg/p^. size of interacting eddy 

relative to particle size x^/k and Reynolds number 

based on the size of the most energetic eddies of 
the flow. An analogous expression to that given by 
Eq . (62) is that based on the Chao’s solution (see 
Chao (1964)). This solution can be considered as a 
substitution of Dingguo’s solution only for k << k 

l.e. for fine particles. 

The last term to be modeled in the VCD group, 
is separated into four correlations: 

18^ 

A 3 = - -J- »S u Li u U + a S u Li u Li ) 
d S 

T) T2 


- (a s u Li u si + a s u u u iii>3 

T3 T4 


(63) 


where the triple correlations T3 and T 4 are modeled 
In a similar manner to that used for the 

calculation of T1 and T2 by using Eqs (23) and 
(24), thus 


p 3Ui‘ .Qq 

T1 = a S u Ll u Li = (k LP /t LP ) (2u Li u Ll “ferV 


12 = “sWl = “Sf>T 


•“Lds, 


T3 = -C 




♦5 (k LP /t LP ){u Li u Si — ftT 


+ Uc.u; 




Si LI dx l J ? 


T4= -C 


^si D s 


du.” , a 


LIS 


>5 (k LT /e LT )(u Li u Sl“5^“ + u Si u Ll J T 

(64) 


These correlations can also be collective 
using single velocity and length-scale, 
k. and t . The first two and the last 


ly modeled 
and total 
two terms 


yield, respectively: 


2 k, 


T1 + T2 " C *5 ^ («, 


A Z 1 

ou. . a 0 
Li S. 

dx, } 


T3 ♦ T4 = C 


2 k. 


1 1 

t t 


♦5 ^ u Li U Si dxj 


(65a ) 
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1.0 


u Si u Li 


- t t 

dx . 


(65b) 


The terms in Group #9 of Eq. A5 are of 
diffusive and dissipative nature. The diffusion 
terms as they appeared in Group #9 are multiplied 
by the molecular viscosity and therefore will be 
neglected due to their relatively small magnitude^ 
Other higher order correlations and mixed 

correlations in this group are also neglected 
according to the modeling principles stated 

previoulsy in section 4.3 and as they are also 
multiplied by the molecular viscosity. 

By substituting all previously modeled terms 
into the exact form of the turbulence kinetic 
energy equation, and rearranging these terms, one 
obtains the simplified modeled form given in 
Appendix B (Eq. B2). 

Since the present model is based on an exact 
equation, namely the turbulence kinetic energy 
equation, and the modeled form of this equation has 
no adjusting coefficients it Is expected that the 
model will generally produce good results and have 
less limitations compared to other models. The 
only modeling assumption is that of the Boussinesq 
gradient type, which generally is accepted. The 
correlations that requires questionable 

semi-empirical modeling assumptions and 

introduction of empirical constant are (i) fewer in 

number (for a^u^u^j and c^Ug^ti^). and (ii) for 
terms having an order of magnitude smaller (by 

ratio a' / cl) compared to other main terms in the 
k-equation. The only significant new correlation 
used in the present closure is that due to the 
relative motion between the phases. This is 
modeled with less restrictions and taking into 
consideration the effect of the particle 
diameter -eddy size ratio on the particle response 
to the eddying motion. The limitation of the 
present one-equation k model closure is the 
algebraic formulation for the length scale. Since 
there are many factors affecting this length scale 
and since it is even difficult in many practical 
applications to give a unique and accurate 
description of the length scale, the use of a 
transport equation for the length scale in the 
two-equation model of Elghobashi and Abou-Arab 
(1983) is expected to give better results with fine 
particles (d g < q). However, it should be noticed 

that the former model and any other similar models 
contain some empirical constants specific for 
various flow conditions, and they require the 
solution for an additional transport equation. It 
can be expected that the present model combined 
with an appropriate length scale equation eg. 
dissipation rate equation will simulate better most 
of the important . features of multiphase turbulent 
flows, particularly the fluid particle interaction. 
In that case the number of closure transport 
equations will increase to three (if two velocity 
scale, kp and Icp, and one length scale transport 

equations) or four transport equations (if two 
velocity scales and two length scale. l p and 

are adopted) . 
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Figure 3. Computed and measured [20] mean velocity 
distribution of air in air-solid pipe 

flow ( — •• single-phase, 2 Eqs. k-t 

model [12], — 1 Eq* k-model [27], -xx- 
1 Eq. k-model [12], 1 Eq . i> t -modei 



“L 
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Figure 4. Predicted radial distribution of the 

turbulence intensity of air in air-solid 
pipe flow using different turbulence 
models (keynote as in Fig. 3). 


7. SAMPLE OF RESULTS 

Figures 3 to 6 compare the present predictions 
with LDA-measurements for single and two-phase 
turbulent, pipe flow (see Maeda et al (1980)) and 
turbulent round water jet laden with uniform-size 
solid particles (see Par thasarathy and Faeth 
(1987)). Both flows are oriented vertically 

downward. These flows are axi syrmetr leal . The 

concentration profiles are given as input data 
based on experimental results. 

In these flow situations the average eddy size 
f was varying from about one dg to few hundreds 
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d g . The corresponding representative eddy size in 

the transfer range was only a fraction of particle 
diameter d g in the pipe flow case. In the Jet flow 

case the mean size of large eddies and the 
Kolmogrov length scale were also varied in a wide 
range. These scales are field variables, and they 
depend upon the flow configuration, location in the 
flow domain and particle dimensions. 

Two-phase flow solutions were obtained by 
solving the flow governing equations in their 
modeled form which are described in the previous 
section and given in Appendix B. The numerical 
procedure used for these predictions is based on a 
developed version of the Gemnix-Code of Spalding 
(1977). However, since the main objective of this 
paper is to give a complete description of a 
developed turbulence closure for multiphase flows 
and due to the space limitation the details of this 
numerical approach will not be given here. The CPU 
Time for the two considered flow cases was about 4 
minutes on a VAX 760 Mini— Compu ter and 2 minutes on 
IBM 3084. 

Case I : Cas-Solid Vertical Pipe Flow: 

Figure 3 shows a comparison between the 
experimental data and the present predictions using 
five different models of turbulence namely ]. 
One-equation k-model of Roco and Hahadevan (1986), 
2. One-equation -model of Roco and Balakrishnan 

(1985). 3. The k-equation as given in the 
two-equation k-e model of Elghobashi and Abou-Arab 
(1983), 4. The two-equation model of Elghobashi and 
Abou-Arab (1983), and 5. The 
present one-equation k-model. The figure displays 
the mean axial velocity distribution In the fully 
developed zone of the pipe flow for single and 
two-phase cases. The differences between the 
predictions of all one-equation turbulence models 
and experiments Is mainly caused by the general 
algebraic expression adopted for the turbulence 
length scale which was not optimized or adjusted. 
In the present computation the concentration 
profiles are assumed based on previous experimental 
data. The inlet concentration distribution is 
taken to be similar to that given by the best curve 
fit after the experimental data of Soo (1967). 
Figure 4 compares the calculated turbulence 

Intensity defined as ju£ 2 (=| l^) / u L with Its 

measured values. The near wall treatment Is based 
on a modified form for the law of wall (see 
Abou-El lai 1 and Abou-Arab (1984), Lee and Chung 
(1987)) and the particle slip condition at the pipe 
wall. 

The present model as compared with the one 
equation models of Refs. [12] and [27] predicts 
slightly higher values for the mean flow 
quantities. However, it must be mentioned that the 
last k-models contain empirical constants In the 
dissipation term of the k equation. Concerning the 
fluctuating flow quantities the present model gives 
slightly better results for the turbulence 
intensity In the near wall region than other 
one-equation models. The expression for the 
turbulence length scale was not optimized. It is 
also expected that the current model will give 
better predictions for coarse (dg > tj) and heavy 



r / z 

Figure 5. Computed and measured [23] mean velocity 
distribution in water-solid jet flow. 



r / 2 


Figure 6. Predicted radial distribution of the 

turbulence Intensity of water in water- 
solid jet flow. 


suspension flows for which no set of comprehensive 
2D data is available for comparison. 

Cs.se II : Turbulent Round Water Jet Laden with 

Uniform-size Solid Particles 

Comparison between experimental data and 
numerical predictions of different turbulence 

closures are given in Figures 5 and 6. 

The two-phase flow measurements on velocity, 
concentration and turbulence correlations used for 
comparisons are taken from Par thasarathy & Faeth 
(1987). Different axial locations within the round 
between eight and fourty Jet diameters from 
the injection nozzle are 
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r/z 


Figure 7. Radial distribution of solid phase 

concentration In water-solid jet flow. 


considered. The radial distribution of solid 
concentration given in Figure 7 are at eight 
diameters from the nozzle. 



0.2KI0' 4 I0* 4 I0* 3 I0' 2 ICf 1 I 10 d s //(5x) 

<b) 


The comparison shows that the two-equation 
two-phase k-t model of Elghobashi and Abou-Arab 
(1983) describes both flows better than the 
one-equation mass/time averaged turbulence model. 
However. It is important to mention here that the 
presently developed closure uses only one transport 
equation and without adjusting any empirical 
coefficient. At the same time, the present closure 
Is in its early stages and more refinements and 
validation tests are required. especially for 
coarse particles two-phase flows for which one 
would expect that the present k-formulation will 
provide improved predictions. The difference 
between the predictions of the mean and fluctuating 
flow velocity components as obtained by the present 
k-formulation and that of Ref. [12] depends on the 
particle size relative to its surrounding eddies. 
Figure 8 illustrates this difference at two 
different loading ratios for the pipe flow at a 
radial distance (R - r)/R equal 0.1. 

6. CONCLUDING REMARKS 

The turbulence closure presented in this paper 
for dilute suspension flow is based on the fluid 
turbulence kinetic energy equation. The main 

features of this model are 

i) Two velocity scales are adopted in computation 
for large and medium size eddies, 
corresponding to the turbulence production and 
transfer range. respectively They are 

expressed into the governing equations by a 
specific local mass/time averaging On this 
basis the spatial and temporal transfer rates 
of the thermodynamic quantities and the 
particle-eddy interaction are better 
estimated. 


Figure 8. Differences between the present 
k-formulation and that of Ref. 

[12] for different mass loading ratios 
L.R. : (a) fluctuating and (b) mean 

axial velocity. 


ii) Spectral analysis of the interaction mechanism 
between particles and most energetic eddies 
provide analytical correlations for closure 
The particle response and the modulation of 
turbulent eddying motion is given as a 
function of the par t ic 1 e- f 1 u id density and 
size ratios. 

lil) To keep the number of transport equations of 
the turbulence closure and the number of 

empirical constants as minimum as possible, 
the length scales ]p and 1^ are described 

using algebraic expressions. Relations 
between these scales (E~;s. 40 and 47) are 

suggested . 

iv) The model does not introduce additional 
empirical constants to the closure of the 
velocity scale equation. 


159 



REFERENCES 


1. Abou-Arab, T.W., "Turbulence Models for 
Two-phase Flows". Encyclopedia of fluid mechanics. 
Vol. 3. Ed. N.P. Cheremi sinof f . Gulf Publ., New 
Jersey (USA), pp. 663-907. 1967. 

2. Abou-Arab, T.W. and Abou-Ellail, M.M.M.. "Heat 
Transfer in Gas-solid Turbulent Pipe Flow". Proc. 
Int . Conf. on Numerical Methods for Transient and 
Coupled Problems. Italy. 19S4. 

3. Abou-Ellail. M.M.M. and Abou-Arab, T.W. , 
"Prediction of Two-phase Flow and Heat Transfer in 
Vertical Pipes". 5th Int. Svmp. Turbulent Shear 
Flows, pp. 8. 1-8.9, 1965. 

4. A1 Taweel . A M. and Landau. J. "Turbulence 
Modulation in Two-phase Jets". Int. Multiphase 
Flow. Vol. 3. pp. 3*11-351. 1977. 

5. Buyevich. Yu. A.. "Statistical Hydrodynamics 

of Disperse Systems". Part 4. Physical background 

and general equations". J. Fluid Mech.. Vol. 40. No 
3. pp. 340-507 . 1971. 

6. Chen, C P, , "Multiple-Scale Turbulence Model 
in Confined Swirling Jet Predictions". AIAA-J . 
Vol. 24. No. 10. 19S6 . 

7. Crowder, R.S.. Doily. J.W., end Humphrey, 

J.A.C.. "Numerical Calculation of Particle 

Dispersion in a Turbulent Mixing Layer Flow", J of 
Pipelines. Vol. 4. No 3. pp 159-170. 1964 

8. Crowe, C.T. and Sharma, M.P. "A Novel 

Physico-computat ional Model for Quasi 

One-dimensional Gas Particle Flows", Trans. ASME. 
J. of Fluids Engineering. Vol. 100, pp. 343-349. 
1978. 

9. Ifanon , H., Wolfshtein. M and Hetsroni, G. 
Numerical Calculations of Two-phase Turbulent 

Round Jets", Int. Multiphase Flow. Vol. 3, pp. 
223-234. 1977. 

10. Dingguo, X.. "Turbulent Kinetic Energy of 

particle Phase in Solid-fluid Suspension 
Turbulence", Int. Symp. Multiphase Flow-China, pp . 
396-401. 19S7. 

11. Elghobashi, S.E and Abou-Arab, T.W . "A 

Second Order Turbulence Model for Two-phase Flows". 
Proc. of the 7th Int. Heat Transfer Conf.. TF4. pp 
219, 1962. 

12. Elghobashi. S.E. and Abou-Arab. T.W . "A 

Two-equation Turbulence Closure for Two-phase 
Flows". Phys. Fluids, Vol. 26. No 4, pp 931-93S 
1963. 

13. Elghobashi. S E. . Abou-Arab. T W . Rizk, M.. 
Mostafa. A., "A Prediction of the Particle Laden 
Jet with a Two Equation Turbulence Model". Int. J 
Multiphase Flow. Vol. 10. No 6. pp 697-710. 1964 

14. Genchiev Ah D and Karpuzov. D.S . "Effects of 

Motion of Dust Particles on Turbulence Transport 
Equations", J. Fluid Mech . Vol. 101. No.4, pp 

633-642. I960. 

15. Hanjalic', K. and Launder. B.E.. "A Reynolds 
Stress Model of Turbulence and its Application to 


Thin Shear Flows. J. Fluid Mech.. Vol. 52. No 4, 
pp. 609-63S (1972). 

16. Hanjalic’, K.. Launder. BE. and Schiestel. 
R.. "Mul t iple-time Scale Concepts in Turbulent 
Transport Modeling”. Proc. of the 3rd Int. Symp . on 
Turbulent Shear Flow, pp . 10-31, 1979. 

17. Hetsroni, G. , "Handbook of Multiphase Flow". 
McGraw Hill. New York. 1982. 

18. Hinze, J.0. . "Turbulence", McGraw-Hill. New 
York. 1975. 

19. Launder, B.E., "Heat and Mass Transport". 
In. Turbulence, Ed. P. Bradshaw, Springler-Ver lag , 
Berlin, pp. 231-287, 1978. 

20. Haeda , M. , Hishida, K. and Furutani, T. , 
"Optical Measurements of Local Gas and Particle 
Velocity in an Upward Flowing Dilute Gas-solids 
Suspension", Proc. Polyphase and Transport 
Techno logy -San Francisco, pp. 211-216, 1980. 

21. Melville, W.K. and Bray, K.N.C.. "A Model of 
the Two-phase Turbulent Jet", Int. J. Heat Mass 
Transfer, 22, pp. 647-656, 1979. 

22. Michaelides, E.E., and Farmer, L.K., "A Model 
for Slurry Flows Based on the Equations of 
Turbulence", ASME-FED Vol. 13, "Liquid-solid Flows 
and Erosion Wear in Industrial Equipment". Ed. M.C. 
Roco, 1963. pp. 27-32. 

23. Parthasarathy , R.N. and Faeth. CM.. 
"Structure of Particle-laden Turbulent Water Jets 
in Still Water", Int. J. Multiphase Flow, 13(5). 
pp. 699-716. 1987. 

24. Peskin, R.L. , "Stochastic Application to 
Turbulent-diffusion". In Int. Symp on Stochastic 
Hydraulics, Ed. C.L. Chiu, Univ. of Pittsburg, PA. 
pp. 251-267, 1971. 

25. Picart, A., Berlemont, A. and Couesbet, G. . 
"Modelling and Predicting Turbulent Fields and the 
Dispersion of Discrete Particles Transported by 
Turbulent Flows", Int. J. Multiphase Flow. 12(2). 
pp. 237-261, 1986. 

26. Roco. M.C. and Bal akr i shnan . N .. 
"Mul t 1 -Di mens ional Flow Analysis of Solid-liquid 
Mixtures", J. of Rheology. Vol. 29. No. 4. pp. 
431-456, 1965. 

27. Roco. M.C. and Mahadevan. S. , "Scale-up 
Technique of Slurry Pipel ines-Par t A' Turbulence 
Modeling". ASME-J . of Energy Resources Technology'. 
Vol. 108. pp. 269-277, 19S6. 

28. Roco, M.C. and Shook C.A., "Modeling Slurry 

Flow: The Effect of Particle Size". Canadian J. 

Chem. Eng., Vol. 61. No. 4. pp. 494-503. 1983. 

29. Roco, M.C. and Shook, C.A., "Turbulent Flow of 
Incompressible Mixtures", J. Fluids Engineering. 
Vol. 107, June 1985. pp. 224-231. 

30. Shuen, J.S. , Chen. L.D. and Faeth. CM, 

"Prediction of the Structure of Turbulent Particle 
Laden In Round Jets", AIAA J., Vol. 21. pp 
14S0-1483, 1983. 


160 



31. Soo. S.L., "Multiphase Fluid Dynamics”, 
Revised Edition, S.L. Soo Assoc., 1983, also "Fluid 
Dynamic of Multiphase Systems". Blaidesll. 
Waltham. MA, 1967. 


- p K^ (u Ki u KJ> t+ VkiJ + (I Ki>-K 

- 7 7 , , Interactions 

Co 1 1 i s i ona Effect with ( „ K) 


(A 2 ) 


32. Spalding. D.B. "A General Computer Program for 
Two-Dimensional Parabolic Phenomena". Dept. of 
Mech. Engng. , Report No. HTS/77/9. Imperial 
College, London, 1977. 

33. Vernier, Ph.. and Delhaye, J.H.. "General 
Two-phase Flow Equations Applied to the 
Thermodynamics of Boiling Water Nuclear Reactors", 
Energie Primaire, Vol. 4. No. 1-2, p. 5. 1968. 


where 


K = phase (or generally a flow component) 

f^ = mass/time average of f over K 

i.J = 1.2,3 (Cartesian coordiantes) 

= body force in the i-th direction 
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APPENDIX - A 


Mass/^Time Averaged conservation Equations 

Mass averaging the conservation equations of mass 
and momentum over- a flow component (K) one obtains 
a new system of equations for mean velocity, 
concentration and kinetic energy of turbulence. 

The point instantaneous conservation equation 
can be written for any flow component (K) or for 
the entire mixture in the folowing general form 

g— (p^/) + v(p^u) + V*J - S = 0 (Al) 


interaction vector 


The interaction term as it stands for 
solid/liquid drag is given by 

. t . , . . _ __ C DS ^ u Ly u Sy ^ u Lv U Sv^ 

^Si^L " ” °' 75a S p L d .. .1.7 

S ' 1 _a S' 


(A3) 

This drag term takes a simple form for Reynolds 
numbers less than unity 

(l S i>L = “ (I L1>S = ^( 18 h/ d sH u u -u sl > <A4) 

The transverse effects caused by the presence of 
other solid particles, Saffman force and Ho and 
Leal inertial force are neglected. Their 

importance is small in dilute suspension turbulent 
flows with fine particles. 


where: p = densi ty 

u^ = velocity vector 

- transported quantity 
2 = flux vector for p 

S = source term 

Let assume * = 11 ^. By splitting each flow 

property into mean and turbulent fluctuating 

component (u^ + u^} and mass/time or double time 

averaging the equation (Al). one obtains the 
following momentum equations in the i-th direction 
for an incompressible phase (K) without mass 
exchange with other flow components (see Roco and 
Shook (1985)): 


P K ft (a K u Ki + a K u Ki ) + P K 3 x (a K u Ki u Kj ) 

- J 

Time Rate Mean Flow Convection 


d “ — “ _ 3 ^ ^ 3 a^ 

^K^K^i ■ ax^ (‘V’K +a K p K ) + Pk aj“ + ^ dx. 

Body Force Pressure Effect 


+ 5 T tVK 2 . f ” p K a K u Ki u Kj 

* - M — 

Frictional Inertial 

Effect Effect 


Equation (A2) contains terms due to the 
unsteady flow, mean flow convection, diffusion, 
pressure, body force, as well as frictional, 
inertial and collisional effects. The mean form of 
the turbulence kinetic energy governing equation is 
obtained by subtracting from the steady state 
ins tantaneous momentum equation for a component K = 
L the corresponding mean flow equations, and then 
multiplying the resulting difference equation by 
(u’ + u" .). By averaging one obtain the kinet ic 

K i K i 

energy equation for a flow component K = L. This 
equation reads 


-- ,^Ki /2 du K? /2 

fWKjt-dXj - + 

Group #1 (Convection) 


> + ^ p kV u Kj u Ki u Ki 

Group #2 (Velocity 

Di f fusion) 
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u,, .u.. . u r 

kj ki ki 


+ u Kj u Ki u Ki + u Kj u Ki u Ki + U Kj U Ki u Ki )] 

+ ( WKi U K; U Kj * Vki u K-. u Kj ) 

Group #3 (Higher Order Correlations) 
+ other 4th order and minor terms)]= 
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" “k ^ {u Ki P K + “ki’k + “ki'k + u "ki P K ) 

Croup #4 (Pressure Diffusion) 
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^ +p ..^ + P^ii + rf!lii 

ax. k ax. k ax< k ax. 

i i i i 


Group #5 (Extra Production and Transfer) 
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- ^“K 9T { * U Kl a K'a^' J + U Kl a K _ STj 


au, au, du” au" 

+ 'YVk ST -at + Vic"* ~aT “AT + mixed 

J J J J correla- 


Croup U9 (Extra Dissipation & Diffusion) 
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Group #10 (Extra Dissipation) 


APPENDIX B 
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au Ki Dux . 
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Croup US (Production) 
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Croup U1 (Extra Production) 
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The M odeled Form of the Turbulence Kinetic 
Energy Equation 


The steady state turbulence kinetic energy equation 
for the liquid phase (K = L) is: 

_ _ ai^ PLt ^li 2 

p L a L u Lj ~ axj ~o^ a~ + VYt ( a^ _) 

Convection Diffusion Production 


lii,!!,- fit, .lift,- 

Vi 8 *! '*>1 *-j dxj> a^Sx, (“ 

Extra Production 


au, 


Li 


dx . ' 
J 


• p l“l £ l 

Di ssipat ion 
+ 18 ^L - - . u Lt 

+ -j- {u Lr u si ) ^r at 
d s ^ 1 

Extra Dissipation 

J8m, 

■ a S k L (, " r SL +r RL ) - — (T1 > T2 - T3 - T4) (Bl) 
d S 

where the expressions for r gL , T1 . T2. T3 and 

T4 are given in the text. 
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